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Advances  in  the  Study  of  Separated  Flows  /I 

Zhuang  Fenggan  (Beijing  Institute  of  Aerodynamics)  and  Zhang 
Hanxin  (Chinese  Aerodynamic  Research  and  Development  Center) 

Abstract 

This  paper  is  divided  into  two  parts.  The  first  part  deals 
with  the  discussion  of  numerical  simulation  of  separated  flow  in 
which  problems  requiring  further  investigation  using  the  NS 
equation  are  identified.  The  second  part  involves  the  use  of 
differential  topology  for  the  quantitative  analysis  of  a  flow 
field.  Furthermore,  some  preliminary  investigation  on  three- 
dimensional  separation  and  vortex  formation  is  discussed. 

I.  Introduction 

For  a  long  time,  people  believed  that  flow  separation  should 
be  avoided  in  order  to  obtain  good  aerodynamic  properties .  In 
reality,  separation  always  occurs.  For  instance,  separation 
always  occurs  on  the  lee  side  when  an  aircraft  flies  at  a  large 
attack  angle.  In  addition,  we  also  discovered  that  the  leading 
edge  separation  vortex _could  improve  the  lift  of  the  airfoil. 

Hence,  the  important  task  is  to  understand  and  control  the 
mechanism  of  separation  and  vortex  motion  in  order  to  create  the 
condition  for  implementing  active  control.  This  paper  will  not 
discuss  any  experimental  techniques.  It  only  discusses  the  most 
important  progress  in  digital  simulation  and  qualitative  analysis 
from  the  theoretical  angle. 


II.  Numerical  Simulation  of  Separated  Flow 

First,  we  must  have  a  mathematical  model.  One  model  is  used 

to  simulate  a  steady  flow  and  the  other  is  to  simulate  a  real 

time  flow.  The  steady  and  non-steady  NS  equation  can  be  used  as 

an  accurate  mathematical  model,  respectively.  This  also  includes 

turbulent  flow.  Of  course,  equations  can  be  simplified  in  many 

cases.  Various  turbulent  flow  models  have  to  be  used  in 

calculating  a  turbulent  flow.  Different  simplications  procedures 

may  be  used  in  different  regions.  In  the  solution  finding 

process,  iterations  can  be  performed  in  order  to  match  it  between 

regions.  This  method  is  ve~y  often  suitable  for  a  steady  flow. 

Even  for  -a  steady  flow,  iterations  between  regions  may  also  be 

lengthy.  Therefore,  we  employed  simplified  NS  equation  which 

satisfy  all  regions  within  a  specific  accuracy  range.  The 

current  simplified  NS  equations  include:  the  simplified  NS 

[  1  ] 

equation  introduced  by  Davis  based  on  Zheng  Xiangi's  viscous 

_  i 

thin  shock  wave  layer  theory,  which  is  accurate'  to  the  Re  order 
of  magnitude  and  is  called  the  viscous  shock  wave  layer  equation, 
the  parabolic  NS  equation,  i.e.,  the  PNS  equation,  introduced  by 
Lubard  and  HelliwellL  J  by  omitting  all  diffusion  terms  related 
to  the  flow  direction  in  the  steady  NS  equation,  and  a  thin  layer 
approximation  equation  which  omits  all  diffusion  terms  in  the 
flow  direction  and  transverse  direction.  In  two-dimensional 
cases,  PNS  equation  is  identical  to  the  thin  layer  equation.  In 
steady  conditions,  the  viscous  shock  wave  layer  equation  is 
parabolic  in  both  the  flow  and  transverse  directions.  However, 
due  to  the  presence  of  a  subsonic  region,  the  equation  is  also 


weakly  elliptical.  PNS  equation  is  parabolic  in  the  flow 
direction.  However,  it  is  elliptical  in  the  transverse 
direction.  It  can  be  used  to  simulate  a  transverse  separation 
problem  where  the  flow  separation  region  is  not  too  large. 
Similarly,  because  of  the  presence  of  a  subsonic  region,  the 
pressure  gradient  along  the  flow  direction  should  be  specially 
treated  in  performing  iterations. 

this  paper  was  received  on  June  26,  1984 

When  a  separation  region  is  present,  a  single  scan  most  often 
cannot  yield  accurate  results.  Therefore,  it  is  necessary  to 
perform  repeated  scans  or  overall  iteration.  A  simplified 
equation  is  only  suited  for  a  separation  region  which  is  not  very 
large.  Normally,  the  complete  NS  equation  will  have  to  be  relied 
on  to  directly  deal  with  the  problem. 

Currently,  the  use  of  NS  equation  (including  simplified  NS 

equations)  to  simulate  a  viscous  separated  flow  is  in  the 

.  ■*. 

developing  stage.  A  great  deal  of  work  has  been  done  on  laminar 
and  turbulent  flows.  Simulation  is  partially  successful  in 
dealing  with  some  typical  flows  such  as  shock  wave  boundary  layer 
interference,  two-dimensional  compression  turning,  flow  around 
the  airfoil  and  flow  around  a  hemispherical  cylinder.  Especially 
in  laminar  flow  conditions,  the  data  are  basically  in  agreement 
with  the  experimental  results.  In  the  case  of  turbulent  flow, 
there  is  usually  a  larger  difference  because  there  is  a  lack  of 
reliable  turbulent  flow  models  in  the  separation  region  and  the 
neighboring  regions.  In  addition,  the  Reynold's  number  used  in 


the  simulation,  which  is  to  some  extent  related  to  the  turbulent 
model,  is  not  high  enough.  Therefore,  current  numerical 
simulation  of  NS  equation  cannot  yet  be  used  extensively  to  solve 
separated  flow  problems  in  engineering  due  to  the  limitation  of 
computer  and  turbulent  model.  In  spite  of  these  limitations, 
numerical  simulation  of  NS  equation  is  very  useful  to  the  further 
understanding  of  the  flow  mechanism.  Furthermore,  it  creates 
conditions  for  engineering  application  on  the  long  run.  In 
steady  conditions,  NS  equation  is  hypobolically  elliptical. 
Because  of  the  elliptical  nature,  its  numerical  simulation  is 
very  complicated.  After  taking  the  development  of  new  computers 
into  consideration,  we  leaned  toward  a  time  correlation  method, 
which  can  reflect  real  time  simulation,  to  solve  NS  equation.  In 
order  to  make  practical  numerical  simulation  feasible,  we  must 
properly  solve  the  following  problems. 

1.  The  coordinate  system  and  mesh  chosen  must  have  a 

sufficiently  high  resolution  to  describe  the  physical 

characteristics  of  the  flow  field.  Currently,  a  body  coordinate 

is  most  frequently  used  to  turn  the  physical  space  to  solve  the 

equation  into  a  simple  regular  region  (such  as  a  rectangular 

parallelepiped  or  a  sphere) .  The  boundary  of  the  physical  space 

r  3  q 

is  the  boundary  for  the  simple  region.  The  early  TTM  method 
using  elliptical  equation  to  solve  the  problem  is  usually 
considered  as  the  most  popular  method.  The  elliptical  nature 
ensures  the  smoothness  of  the  solution.  Furthermore,  we  can 
adjust  the  distribution  of  the  nodal  points  of  the  mesh  by 


choosing  an  appropriate  control  function.  The  disadvantage  of 

this  method  is  that  it  takes  more  computer  time  and  memory  space. 

A  simpler  approach  is  to  use  a  mesh  formation  method  for 

[4] 

algebraic  equations  .  A  typical  algebraic  method  is  to 
introduce  several  layers  of  control  surfaces  between  the  inner 
and  outer  boundary  based  on  the  characteristics  of  the  physical 
problem  so  that  the  boundary  corresponds  to  plane  in  the  simple 
region.  Then,  the  corresponding  relation  of  internal  nodal 
points  is  found  by  intrapolation  of  algebraic  expressions.  As 
for  recent  advances,  in  addition  to  using  various  methods  to 
generate  meshes,  the  study  of  optimization  of  meshes  and  self- 
adapting  .meshes  is  more  important  which  allows  the  automatic 
capture  of  shock  wave  and  vortex  surface.  It  can  accurately 
reflect  the  spatial  flow  in  the  large  gradient  region.  The  key 
is  to  correctly  present  a  guideline  for  optimization  and  self¬ 
adaption.  This  is  a  problem  worth  further  investigation, 

2.  The  selection  of  a  difference  scheme  can  directly  affect 
the  stability  and  rate  of  convergence  of  the  solution.  Explicit 
time  dependent  methods  were  used  in  earlier  simulation.  However, 
due  to  limitation  of  stability  conditions,  the  allowed  time  step 
in  the  calculation  is  too  small.  Therefore,  many  implicit  methods 
have  been  developed.  Here,  only  implicit  difference  schemes  are 
analyzed . 

First,  the  center  of  our  discussion  is  implicit  difference 
methods,  which  include  the  non-iterative  implicit  factorization 
method^^^,  ADI  method  and  time  sharing  implicit  center 
difference  methodL  .  It  was  proven  that  these  three  methods  are 


interchangeable  within  the  range  of  second  order  precision. 
Moreover,  all  three  use  center  difference.  These  are  feasible 
methods.  The  Courant  number,  which  is  related  to  the  time  step 
length,  can  reach  an  order  of  magnitude  of  several  hundreds  or 
thousands.  However,  in  the  calculation,  especially  in  dealing 
with  problems  with  hidden  shock  wave,  it  is  necessary  to  add  an 
artificial  dissipation  term.  This  is  a  key  technical  measure  to 
control  oscillation  and  it  is  highly  empirical.  In  addition,  it 
is  necessary  to  form  viscous  and  inviscous  Jacobian  matrices  and 
to  solve  three  diagonal  matrices  by  chasing  in  the  calculation. 
Hence,  the  computation  and  programming  workload  is  large. 

.  Another  important  implicit  scheme  was  introduced  by 
rgi 

MacCormack  in  1981  .  As  compared  to  the  method  described 

above,  it  has  the  following  advantages.  First,  the  difference 
scheme  employed  is  a  two-point  single  side  difference. 

Therefore,  it  is  only  necessary  to  chase  two  diagonal  matrices. 
Second,  it  is  only  necessary  to  calculate  the  inviscous  Jacobian 
matrix  in  the  computation.  However,  the  effect  of  viscosity  is 
taken  into  account  in  the  stability  condition  when  we  calculate 
the  absolute  matrix.  Third,  when  the  time  step  satisfies  the 
explicit  expression  stability  condition,  it  is  able  to 
automatically  change  into  the  early  stage  MacCormack  explicit 
scheme.  These  special  features  reduce  the  workload.  However,  we 
must  still  calculate  the  characteristic  value  of  the  Jacobian 
matrix.  Moreover,  it  is  necessary  to  create  an  absolute  matrix. 
This  cancels  some  of  the  advantages.  Unfortunately,  an  empirical 
artificial  dissipation  term  must  still  be  added  to  eliminate 


6 


oscillation  in  the  calculation  of  a  separated  flow  field  with 
shock  waves . 

In  order  to  solve  this  problem,  Coakley  ^introduced  an 
implicit  windward  difference  method.  As  a  matter  of  fact  this 
method  originates  from  the  separation  coefficient  matrix  method 
used  to  solve  hypobolic  equations.  The  basic  idea  is  to  use  the 
characteristic  value  of  the  inviscous  portion  of  the  Jacobian 
matrix  in  the  NS  equation  to  control  the  spatial  derivative  of 
the  inviscous  part.  When  the  characteristic  value  is  greater 
than  zero,  a  second  order  difference  directed  toward  the  rear  is 
used.  When  the  characteristic  value  is  less  than  zero,  a  second 
order  difference  directed  forward  is  used.  The  spatial 
derivative  of  the  viscous  portion  employs  center  difference. 

This  scheme,  in  theory,  is  dissipating.  There  is  no  need  to  add 
another  dissipation  term.  Examples  of  shock  wave  and  boundary 
layer  interference  calculation  show  that  this  method  is  highly 
accurate  and  stable.  The  time  step  for  calculation  is  large. 

The  convergence  rate  is  fast.  Furthermore,  it  does  not  need  an 
added  artificial  dissipation  term.  A  detailed  analysis  indicates 
that  the  region  where  the  characteristic  value  changes  sign 
indicates  that  the  region  where  the  characteristic  value  changes 
sign  may  require  special  treatment.  Coakley  was  inspired  by 
Harten's  work  to  create  dissipation  functions.  These  are 

problems  to  be  studied  further.  Other  difference  schemes  are  not 
discussed  individually. 

3.  In  implicit  calculation,  it  is  very  important  to  provide 
the  computation  condition  on  the  boundary.  If  it  is  not  properly 


given,  the  computation  is  not  going  to  be  stable.  Or,  this 

difference  equation  set  and  boundary  condition  make  the  problem 

unsuited  to  define.  As  an  example,  we  are  quoting  the  paper  by 

[12] 

Abarbanel  and  Murman  .  They  gave  some  inspiring  conclusions. 
In  the  following  equation 

dU  _  dFOJ )  dGOJ) 

dt  dx  dy 

if  the  region  for  calculation  is  0^x^®,-«»<y<<»,  t>0,  they  proved 
that  if  the  backward  Euler  implicit  scheme  or  Crank-Nicolson 
implicit  scheme  is  used  to  solve  the  above  equation,  then 

(1)  the  computation  is  stable  if  the  boundary  condition  is 

Utl'mtUW -u\y 

(2)  the  computation  is  stable  if  the  boundary  condition  is 

Uril=2U\,-U\;' 

Here,  the  first  subscript  represents  the  nodal  point  in  the  x- 
direction  where  "0"  represents  x  =  0.  The  second  subscript  is 
the  nodal  point  number  in  the  y  direction.  It  is  easy  to  see 
that  the  first  boundary  condition  is  implicit  and  the  second 
boundary  condition  is  triply  explicit. 

4.  Accelerated  convergence  in  the  computation  is  also  an 
important  problem.  There  are  many  measures  introduced  in  the 
literature.  We  believe  that  we  should  pay  more  attention  to 
multi-mesh  and  the  self-adapting  mesh  technique  mentioned  above. 

The  multi-mesh  technique  was  first  introduced  by 
in  1962^^.  Later,  Brandt^  ^  made  further  advances.  Its  basic 


advantage  is  that  the  coarse  mesh  can  eliminate  errors  caused  by 
low  wave  number  and  the  fine  mesh  can  eliminate  errors  caused  by 
high  wave  number.  This  technique  has  already  successfully  been 
used  in  the  calculation  of  the  transonic  small  perturbation  equation, 
complete  velocity  potential  equation  and  Euler  equation.  In 
addition  it  is  beginning  to  be  applied  to  the  computation  of  NS 
equation.  Stubbs^ ^  combined  this  technique  with  MacCormack's  /4 
implicit  scheme  and  obtained  very  good  results.  In  the 
following,  we  will  use  a  one-dimensional  problem  as  an  example  to 
split  the  NS  equation  into  inviscous  and  viscous  portions  and  to 
employ  the  method  suggested  by  Stubbs  to  explain  the  procedure  of 
multi-mes.h  calculation  for  the  inviscous  portion.  Let  us  assume 
that  the  inviscous  part  of  the  equation  is 

dU  dF  r 

-dr^-dr~G 


The  procedure  is : 

(1)  to  use  MacCormack's  implicit  method  to  calculate  6UV  = 
(Uj+1-Uj)h  on  the  fine  mesh  nodal  points  (assuming  the  mesh 
spacing  is  h). 

(2)  to  use  the  Ni  scheme  to  calculate  6Uf  on  the  nodal 
points  of  the  coarse  mesh  (assuming  the  mesh  spacing  is  2h  and 
nodal  points  are  j,  j+2-,  j+4,...) 

where  A  =  aF/aU,  B  =  aG/aU.  According  to  this  result,  the  values 
on  other  nodal  points  of  the  original  fine  mesh  may  be  obtained 
by  linear  intrapolation . 


(3)  If  the  mesh  is  still  not  coarse  enough,  we  can  use  the 
Ni  scheme  to  calculate  the  values  on  nodal  points  of  a  coarser 


mesh  (j,  j+4,  j±8,...): 

2 

Values  on  nodal  points  of  the  fine  mesh  can  be  obtained  by  linear 
intrapolation. 

(4)  to  use  MacCormack's  method  to  perform  implicit  fine  mesh 
calculation  based  on  the  values  obtained  with  the  coarse  mesh 
nodal  points.  Then,  steps  (2)  and  (3)  are  repeated  until  a 
stable  solution  is  obtained. 

Of  course,  the  problems  discussed  above  are  the  major  ones 
to  be  seriously  studied  if  numerical  simulation  NS  equation  of 
separated  flow  is  going  to  become  practical.  To  some  extent, 
these  are  fundamental  problems.  On  one  hand,  we  need  theoretical 
studies.  On  the  other  hand,  we  must  rely  on  a  great  deal  of 
numerical  experiments  to  solve  these  problems.  As  a  part  of  the 
theoretical  work,  differential  topology  will  be  applied  to  the 
qualitative  analysis  of  the  flow  field. 


III.  Qualitative  Analysis  of  the  Flow  Field 
Separation  starts  from  the  surface  of  an  object.  Then,  the 
vortex  penetrates  deep  -into  the  flow  region.  The  study  of 
separation  can  begin  with  the  study  of  the  flow  pattern  on  the 


surface  of  the  object.  Lighthill 


[16] 


is  the  pioneer  in  this 


field.  He  introduced  the  friction  line  concept  and  assumed  that 
friction  is  a  continuous  two-dimensional  vector  field.  We 
usually  believe  that  surface  friction  line  is  the  limiting  stream 
line.  If  an  x,y  coordinate  is  introduced  to  the  surface  and 


t  and  o  are  the  friction  components  in  x  and  y  direction, 
respectively,  then  the  differential  equations  defining  the 
friction  line  are: 

Or.tfi  Re,M,a,— ) 

°  (x,yt  Re,  M,a,  ••• ) 

where  Re,  M,  a...  are  parameters  such  as  Reynolds  number,  Mach 
number  and  attack  angle,  and  X  is  a  chosen  variable.  Relative  to 
different  parameters,  the  distribution  of  the  surface  friction 
line  is  different,  i.e.,  with  various  types  of  topologic 
structures.  From  the  qualitative  theory  of  differential  equation 
we  know  that  the  topological  structure  is  determined  by  the 
number  of  singular  points  (points  where  t=o=0),  the  property  of 
each  singular  point  and  the  connection  of  mutual  friction  lines 
between  singular  points.  Furthermore,  according  to  qualitative 
theory,  a  structurally  stable  singular  point  can  only  be  a  first 
order  singular  point.  One  type  is  the  settle  p'oint  whose 
Poincare  index  is  -1.  The  other  type  includes  focal  point,  nodal 
point  and  center  point,  and  the  corresponding  Poincare  index  is 
+1.  For  an  enclosed  curve  surface,  the  total  singular  point 
index  S  is 

S  =  2-2g 

where  g  is  the  topological  seed  value  of  the  curve  surface.  For 
a  curve  surface  topologically  equivalent  to  a  spherical  surface, 
g=0,  g=1  if  it  is  equivalent  to  an  annular  curve  surface. 
Therefore,  there  are  2  more  nodal  points  (including  focal  points 
and  centers)  than  settle  points  on  a  spherical  surface.  Hence, 


their  topological  structure  can  only  be  N2,  N^,  ,  N^,  S£,  N^, 

. ^n^n-2*  w^ere  N  represents  a  nodal  point  and  S  represents 

settle  point.  The  subscript  represents  the  number  of  that 
particular  point.  Bippes^^  obtained  a  dozen  or  so  different 
flow  patterns  up  to  in  the  hemispherical  cylinder  oil  flow 

experiment.  We  are  not  going  to  investigate  very  complicated 
flows  here.  Instead,  we  will  use  the  flow  around  an  ellipsoid 
with  an  attack  angle  studied  by  Wang  Guozhang  as  a  typical 
example  for  analysis.  In  addition,  we  hope  to  clarify  some 
controversies.  Based  on  experimental  results  and  theoretical 
analysis,  Wang  divided  the  basic  flow  into  the  following  three 
regions: 

Cl)  0<a<agv,  closed  separation  emerges  at  the  tail  as  shown 
in  Figure  (la) ; 

(2)  agv<a<aav,  open  separation  begins  on  the  leeward  side, 
together  with  the  formation  of  a  small  second  order  closed 
separation  region  in  the  rear,  as  shown  in  Figure  (ifc); 

(3)  a>aaV,  another  closed  separation  appears  on  the  leeward 
side  as  shown  in  Figure  (1c). 
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Figure  1.  Basic  Patterns  of  Flow  Around  an  Ellipsoid  with  an 

Attack  Angle 

1.  separation  line 

2.  side  view 

3.  top  view 

4.  open  separation  line 

5.  second  order  separation  line 

6.  side  view 

7.  top  view 

8.  separation  line 

9.  side  view 

10.  top  view 

The  boundary  of  each  region  shown  in  the  figure  is  related 
to  the  aspect  ratio  of  the  ellipse  when  the  incoming  flow  is 
identical.  This  flow  pattern  is  quite  different  from  the  flow 
around  a  slender  body  at  a  large  attack  angle.  Here,  the  shape 
of  the  head  and  the  tail,  especially  the  head,  has  an  important 
effect.  Before  analyzing  these  patterns,  several  points  are 
obvious  from  the  topological  viewpoint.  The  separation  line  and 
re-adhesion  line  must  be  the  limiting  stream  lines.  The  limiting 
stream  lines  on  both  sides  of  the  separation  line  will  be  drawn 
closer  to  the  separation  line.  The  limiting  stream  lines  on  both 


sides  of  the  re-adhesion  line,  however,  are  moved  away  from  it. 
The  separation  line  may  be  a  part  of  the  limiting  stream  line  or 
may  be  the  entire  limiting  stream  line.  A  limiting ' stream  line 
must  be  from  a  singular  point  to  a  singular  point.  Or,  it  must 
be  an  enclosed  curve.  It  cannot  terminate  at  a  regular  point. 


With  regard  to  closed  separation,  there  is  no  controversy. 

Figure  (la)  shows  that  the  separation  line  is  the  line  of 

friction  from  S  to  Ng ,  which  agrees  with  the  N^S-j  topological 

structure  model.  The  open  separation  Wang  introduced  to  Figure 

(it)  was  in  fact  first  observed  experimentally  by  Maskell. 

Wang  pointed  out  that  this  open  separation  is  a  separation 

without  any  singular  point.  It  is  not  related  to  Lighthill's 

singular  point  separation  at  all.  In  addition,  Wang  also  made 

a  transition  from  closed  separation  to  open  separation  and 

r  ion 

attributed  it  as  the  so-called  tongue  splitting  .  These 
conclusions  warrant  further  investigation  because  tongue 
splitting  does  not  agree  with  laws  of  topology..  In  addition,  as 
shown  in  Figure  (1b),  the  open  separation  line  must  be  located  on 
the  line  of  friction  connecting  ,  to  Ng.  If  the  line 


/6 


terminates  at  a  new  N  ,  then  it  obviously  does  not  satisfy 

9 


Poincare's  law  of  topology.  Open  separation  is  an  important 
three-dimensional  separation  mode  from  the  formation  of  vortex. 
There  is  practical  significance  to  study  it  further.  We  noticed 
that  Tai^^l  pointed  out  that  drawing  limiting  stream  lines 
together  alone  is  not  enough  to  create  vortex  separation.  It 
also  requires  "collision"  with  the  inviscous  stream  line  on  the 
outer  fringe  of  the  boundary,  just  as  the  pattern  shown  by 
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Maskell  earlier.  Thus,  we  analyze  that  separation  and  vortex 

[21  ] 

formation  must  be  studied  as  a  whole.  Dallmann  first  pointed 

out  that  the  spatial  topological  structure  corresponding  to  an 

identical  surface  flow  pattern  is  not  unique.  We  still  do  not 

have  a  complete  topological  theory  in  a  three-dimensional  space. 

However,  if  we  consider  any  surface  in  the  space  and  the 

continuous  velocity  field  on  the  surface,  this  surface  can 

intersect  with  the  surface  of  the  object.  By  defining  the 

singular  point  on  the  body  surface  is  a  semi-nodal  point  N'  or 

[221 

semi-settle  point  S',  then  Hunt  proved  that 


«-( 

where  ZN  is  the  total  nodal  point  number,  z<,  is  the  total  settle 
point  number,  Z^tis  the  total  semi-nodal  point  number,  Zg,is  the 
total  semi-settle  point  number  and  n  is  the  chosen  degree  of 
surface  connection.  For  a  single  connected  region  n=1.  For  a 
double  connected  region,  n=2. 

This  law  can  be  used  to  preliminarily  analyze  the  basic  flow 


pattern  around  a  slender  body  at  a  large  attack  angle.  Nielsen 


[23] 


and  Ericsson 


[24] 


summarized  the  experimental  results  from  low 


velocity  to  a  transonic- flow  with  a  transverse  Mach  number  Mc  of 


less  than  0.5.  According  to  the  size  of  the  attack  angle  t, 
the  state  of  flow  can  be  divided  into  five  regions: 


(1)  0<a<a  ,  the  flow  is  adhered  to  the  body  as  shown  in 

S  V 


Figure  (2a); 

<2)  asv 

leeward  area  as  shown  in  Figure  (2b); 


(2)  agv<a<aav, symmetric  stationary  vortices  appear  in  the 
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(3)  aav<a<a^v, asymmetric  stationary  vortices  appear  in  the 
leeward  area  as  shown  in  Figure  (2fl>); 

(4)  a^v<a<a^v, alternating  non-s tationary  vortices  appear  in 
the  leeward  area  as  shown  in  Figure  (2d);  and 

(5)  a>auv,  the  leeward  area  becomes  non-stationary  vortices 
and  tailing  of  random  vortex  motion. 
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Figure  2.  Basic  States  of  Flow  Around  a  Slender  Body  at  Large 
Attack  Angle  (in  A-A  cross-section) 

The  flow  at  M  <0.5  is  often  called  a  subcritical  flow, 
c 

Conclusive  opinions  on  supercritical  flow  at  Mc  ^0.5  are  still  to 
be  further  studied.  The  classification  into  five  regions  is  also 
a  rough  one.  The  actual  flow  is  far  more  complex  than  those 
drawn  here.  A  topological  structure  analysis  of  the  flow  pattern 
shows  that  the  transition  from  one  region  to  another  is  an 
important  subject.  First,  let  us  look  at  Figure  (2a).  There  are 
only  two  semi-settle  points  on  the  plane,  G=-1,  which  agrees  with 
the  topological  law.  Figure  (2b)  only  shows  4  semi-settle  points 
and  2  nodal  points.  Obviously,  there  is  another  settle  point  in 
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the  upper  flow  field.  The  topological  format  is  Figure 

(2c)  shows  the  generation  of  an  asymmetric  lateral  force  on  the 

slender  body  in  region  (3).  Before  it  is  developed  into  an 

asymmetric  vortex,  a  second  order  separation  and  interaction 

between  the  primary  separation  vortex  and  the  secondary  vortex 

exist  on  the  body  surface.  In  the  early  stage,  some  people  / ' 

believed  that  asymmetric  vortices  are  created  by  asymmetric 

separation  on  the  leeward  side.  However,  G.  Chapman  et  al  found 

experimentally  in  1975  that  asymmetric  vortices  are  not  primarily 

created  by  asymmetric  separation.  Instead,  they  are  due  to  fluid 
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dynamic  instability.  Peake  and  Tabak  believed  that  the 
asymmmetry  of  a  fluid  is  related  to  the  instability  of  the 
velocity  cross-section  at  the  upper  settle  point  on  the  cross- 
section  of  the  object.  Furthermore,  this  mechanism  was 
qualitatively  interpreted  based  on  its  topological  structure. 

The  study  on  the  control  and  breakup  of  asymmetric  vortex  is 
still  in  its  infancy. 

Now  let  us  give  a  preliminary  qualitative  analysis  on  the 
spatial  flow  field.  The  Cartesian  coordinate  Xj,  j  *  1,  2,  3  is 
used.  Uj  is  the  velocity  component  in  the  Xj  direction.  Let  a 
singular  point  be  the  origin.  Its  neighborhood  is 


,,_0”  *' »  a”~(d^,) 


The  stream  line  equation  can  be  written  as 


=  AX 


where 


Am  Om#*  X  =*  (Jf, ,xltX\) r„ 


Let  us  assume  the  eigen  value  of  A  is  Xj ,  j  =  1,2,3  and  the 

eigen  vector  corresponding  to  X^  is  .  In  this  case,  the  eigen 
T 

value  of  A  is  also  A  ^  .  The  corresonding  eigen  vector  is  B^. 
Furthermore , 

BijCkj  =6ik 

The  solution  of  the  differential  equation  is 


f.  =  2  V'*-"  **’’ 


<p,  i  =  1,2,3  are  integration  constants.  X  satisfies  the 
following  third  order  equation 


A~[d(xltxa)  d(xirx,)  d(x„x,)  }  dx. 


a:-A:  =  0 


where  A  is  the  value  of  the  determinant  corresponding  to  A.  For 
steady  flow  and  incompressible  flow,  because  of  the  continuity 


equation . 


2 


When  the  three  roots  are  all  real,  they  cannot  simultaneously  be 

positive  or  negative.  This  means  that  singular  points  cannot  all 

be  nodal  points  or  settle  on  three  mutually  perpendicular  planes. 

Furthermore,  we  can  define  three  plane 

B. .  x.  =0 

xj  J 

These  three  planes  are  tangent  to  regular  flow  surfaces.  Where 
one  root  is  real  and  the  remaining  two  are  complex  conjugate 
roots,  it  is  only  possible  to  define  one  plane  tangent  to  the 
regular  flow  surface.  In  this  case,  there  is  a  singular  stream 
line.  Other  stream  lines  are  cylindrical  spinals  around  this 
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singular  stream  line.  Using  the  'Candan  equation  and  assuming 

d<*„*,)  d(x,,x,)  d(xt,x,) 

then  the  necessary  and  sufficient  condition  for  a  pair  of 
conjugate  solutions  is 

The  ordinary  vortex  we  are  referring  to  (including  the  vortex 
created  by  an  open  separation)  is  this  type  of  spiral  flow.  If 
our  only  objective  is  to  determine  whether  there  is  a  vortex 
structure  around  a  specific  point,  we  may  choose  a  coordinate 
system  which  is  moving  with  the  velocity  of  that  point  and  the 
above  analysis  is  still  valid.  By  using  the  results  discussed 
above,  it  is  very  easy  to  prove  that  the  condition  to  have  a 
vortex  leaving  a  symmetric  flow  plane  is  that  the  initial 
singular  point  should  be  a  focal  point  and  the  core  stream  line 
of  the  vortex  (singular  stream  line)  is  perpendicular  to  the 
symmetry  plane.  A  similar  method  can  be  used  to  study  the 
spatial  flow  near  the  surface  which  demonstrates  that  the 
formation  of  a  horn  vortex  originates  from  a  focal  point  singular 
point  on  the  surface.  Incidentally,  we  noticed  that  the  A  matrix 
is  a  real  symmetric  matrix  for  an  imitational  motion.  All  its 
eigen  values  are  real.  There  is,  of  course,  no  such  vortex 
structure . 
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IV.  Conclusions 


This  paper  does  not  involve  any  practical  engineering 
problems  with  complicated  separation  and  vortex  motion.  The 
difficulties  in  theoretical  calculation  are  apparent.  To  date, 
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Graham  and  Hankey'sL  J  calculation  of  asymmetric  vortex  is  most 

[27] 

successful.  However,  from  Lamont's  experiment  ,  the  magnitude 
of  the  lateral  force  is  also  related  to  the  tumbling  angle  and 
Reynolds  number.  Moreover,  in  turning  point  situation,  the 
lateral  force  is  almost  decreased  to  zero.  The  contribution  of 
laminar  flow  separation  to  the  magnitude  of  lateral  force  is 
comparable  to  that  of  turbulent  flow  separation.  However,  the 
distribution  is  not  consistent.  In  addition  to  carrying  out 
detailed  experimental  studies  on  the  rolling  of  spatial  vortex, 
the  mutual  interaction  among  vortices  from  one  topological 
structure  to  another,  and  the  appearance  of  a  random  tail  flow, 
theoretical  qualitative  analysis  will  provide  us  with  great 
assistance.  This  discussion  covered  in  this  paper  only 
classifies  the  status  and  basis  for  further  studies. 

In  the  preparation  of  this  paper,  Comrade  Li  Suxun  was  also 
very  helpful.  The  authors  wish  to  express  their  gratitude. 
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ADVANCES  IN  THE  STUDY  OP  SEPARATED  FLOWS 
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Abstract 

This  paper  consists  of  two  parts.  The  first  deals  with  the  numerical 
simulation  of  separated  flows  and  problems  need  further  study  are  indi¬ 
cated  with  NS  equations  model.  The  second  involves  an  application  of 
differential  topology  to  qualitative  analysis  of  flow  fields  and  ibe  conne¬ 
ction  between  3D-separation  and  vortex  production  is  discussed. 


Some  Problems  in  Discrete  Vortex  Numerical  Modelling  of 
Vortex  Motion  Behind  a  Circular  Cylinder 
Ling  Guocan  (Institute  of  Mechanics,  Academia  Sinica) 

Abstract 

This  is  a  brief  review  of  the  four  problems  encountered  in 
the  discrete  numerical  modelling  of  unsteady  flow  and  vortex 
motion  behind  a  circular  cylinder,  i.e.,  the  discrete  method  and 
the  stability  of  the  vortex  sheet,  boundary  layer  separation  and 
the  positions  of  nascent  vortices,  the  secondary  vortex  problem 
and  the  reduction  of  circulation  of  vortices  in  the  wake.  The 
paper  also  presents  some  topics  which  require  further  study. 

Introduction 

For  an  abruptly  started  unsteady  constant  velocity 
cylindrical  motion  (corresponding  to  an  unsteady  uniform  flow 
around  a  cylinder).,  the  numerical  modelling  under  subcritical 
conditions  at  high  Reynolds  numbers  needs  another  theoretical 
approximation  model  because  the  results  of  the  numerical  solution 
of  the  existing  N-S  equation  are  not  dependable.  The  important 
characteristics  of  the  viscous  separation  flow  of  a  non-stream 
line  convex  object  are  the  formation,  concentration,  offset,  and 
consumption  of  a  vortex.  In  a  two-dimensional  incompressible  N-S 
equation  (expressed  by  £and  <(»),  the  distribution  of  the  vortex 
function  £  with  time  and  space  will  determine  the  entire  flow(<|<). 
A  simplified  theoretical  model  is  based  on  how  to  relatively 
accurately  reflect  the  distribution  of  £.  For  a  high  Reynolds 


number  flow  with  a  valid  boundary  layer,  this  type  of  small 
viscosity  fluid  flow  can  be  expressed  by  using  the  inviscous  flow 
theory  together  with  the  assistance  of  boundary  layer  flow, 
circulation  and  concentrated  vortex.  In  the  flow  field, 
continuously  or  concentrated  vortices  are  expressed,  in 
approximation,  by  many  discrete  small  vortices  of  some  strength. 
This  is  the  discrete  vortex  method.  In  recent  years,  there  are 
significant  advances  in  the  discrete  vortex  method  in  conjunction 
with  the  boundary  layer  theory.  Some  numerical  experiments  on 
the  circular  cylinder  flow  were  successful  in  calculating  the 
disengagement  frequency.  This  method  is  also  very  simple.  The 
dimensionless  time  in  the  conputation  may  be  as  long  as  over  200, 
which  is  equivalent  to  travelling  a  distance  more  than  200  times 
the  radius  of  the  circular  cylinder  after  the  motion  begins. 
However,  because  this  flow  is  relatively  complex,  the  changes  of 
flow  with  time  include  unsteady  separation  of  the  boundary  layer, 
formation  and  variation  of  the  shear  sheet,  creation  and 
development  of  vortices,  alternating  disengagement  of  vortices 
and  the  variation  of  vortices.  To  date,  the  understanding  of  the 
mechanisms  and  patterns  of  some  problems  such  as  the 
characteristics  of  the  rear  shear  sheet  layer,  the  disengagement 
process  of  vortices  and  the  dissipation  of  vortices  in  the  wake 
is  still  inadequate.  There  are  also  some  problems  in  the 
numerical  modelling  using  the  discrete  vortex  method,  which 
require  further  discussion  and  resolution.  The  following  is  a 
brief  review  of  these  four  problems. 
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I.  Discrete  Method  and  Stability  of  Vortex  Sheet  Motion 
The  motion  of  the  vortex  sheet  is  mathematically  an  initial 
value  problem.  However,  the  question  of  the  existence  of  the 
solution  and  its  uniqueness  remains  unproven. 


This  paper  was  received  on  June  14,  1984. 

This  paper  was  a  part  of  the  National  Separation  Vortex  Flow 
Academic  Report  held  on  April  10-20,  1984. 


In  the  actual  computation  of  the  discrete  vortex  method,  the  /I 

point  vortex  group  tends  to  move  randomly  when  the  computation 
gets  too  long.  In  Reference  [1],  the  cause  of  this  random  motion 
was  analyzed  in  the  domain  of  inviscous  flow  theory.  The  induced 
conjugate-  complex  velocity  q(z^)  at  any  point  z^  in  space  by  a 
continuous  vortex  sheet  with  an  intensity  y(.S)  should  be 


Y(St)dSt 


When  z^  is  located  on  the  vortex  sheet,  singularity  is  shown  in 
the  integration.  -The  principal  Cauchy  value  should  be  used.  If 
the  vortex  sheet  is  divided  into  many  small  segments  and  the 
nodal  point  in  each  segment  is  expressed  by  k+(1/2),  then  the 
complex  velocity  can  be  approximately  written  as 


In  |  +$)  '  M 


The  induced  velocity  by  N  discrete  point  vortices  at  point  Zj  is 


?  <*/> 


By  comparing  the  above  we  know  that  the  velocity  of  a  point 
vortex  calculated  by  the  discrete  vortex  method  is  different  from 


the  velocity  induced  by  the  actual  continuous  vortex  sheet.  The 

difference  is  a  logarithmic  term.  Both  are  in  agreement  when 

point  Zj  is  located  at  the  center  of  the  vortex  segment.  The 

logarithmic  term  gets  larger  when  the  distances  between  z^  and 

the  two  ends  are  unequal.  The  discrete  method  was  used  in 

references  Cl]  and  C 2 ]  to  rearrange  the  positions  of  discrete 

conservation  principle  to  make  the  new  position  of  each  vortex  to 

be  at  the  center  of  the  vortex  segment.  Then,  the  motion  of 

these  vortices  was  calculated.  Thus,  the  randomness  of  the 

motion  of  vortices  over  a  long  period  of  time  is  significantly 

decreased.  The  above  discussion  is  based  on  the  inviscous  flow 

theory.  -On  the  other  hand,  when  two  vortices  are  approaching 

each  other  and  the  spacing  is  6<<1,  the  order  of  magnitude  of  the 

induced  velocity  given  by  the  inviscous  flow  theory  is  (1/6)-®. 

In  order  to  eliminate  the  possibility  of  introducing  a  vortex 

nucleus  to  the  singularity,  the  induced  velocity  in  the  vortex 

nucleus  range  can  be  calculated  based  on  the  viscous  flow 
T31 

theory  .  In  order  to  facilitate  engineering  calculations, 
induced  velocity  in  the  vortex  nucleus  can  be  calculated  based  on 
the  Rankine  vortex  model.  The  velocity  can  even  be  made  to  be 
zero^^.  The  use  of  al'l  above  methods  can  reduce  the  random 
motion  of  point  vortices. 

Stability  analysis  of  the  motion  of  a  non-homogeneous  vortex 
sheet  in  contact  with  the  boundary  layer  is  more  complex.  As  a 
vortex  spiral  develops  towards  downstream,  the  intensity 
decreases  with  the  expansion  of  the  vortex  sheet.  Therefore,  it 
is  somewhat  stable  with  respect  to  Helmholtz  stability.  With 
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respect  to  the  Tollmien-Schlichting  instability  of  the  portion  of 
the  shear  sheet  where  the  vortex  is  in  contact  with  the  boundary 
layer,  the  analysis  is  still  inadequate.  This  type  of 
instability  must  be  predicted  in  the  vortex  developing 
process^^ . 

II.  Boundary  Layer  Separation  and  Position  of  Nascent  Discrete 
Vortices 

The  flow  around  a  circular  cylinder  has  a  very  short  initial 
unsteady  period.  Because  there  are  two  different  separation 
guidelines,  i.e.,  Prandtl  and  M.R.S.  guidelines,  the  position  of 
the -boundary  layer  is  different.  The  initial  circumferential 
separation  differs  by  40°  and  the  initial  dimensionless 
separation  time  (t  =  Vt/R)  differs  by  0.3.  Correspondingly,  the 
position  and  time  of  the  initial  separation  vortex  are  different. 
Reference  [4]  pointed  out  that  it  is  a  reasonable  consideration 
to  place  the  position  of  a  newly  created  vortex"  at  the  location 
where  the  boundary  layer  begins  to  thicken  significantly  and 
separation  bubble  begins  to  emerge  based  on  the  M.R.S.  guideline. 
This  paper  also  compared  the  effect  of  both  situations  on  the 
initial  wake.  After  a  very  short  period  of  time,  when  *t  -1 ,  the 
variation  of  the  cylindrical  boundary  layer  tends  to  vary  in  a 
quasi-steady  state.  Afterwards,  there  is  no  obvious  difference 
between  the  boundary  layer  separation  positions  predicted  by  / 1 2 

these  two  guidelines.  Thus,  the  early  stage  of  variation  at 
t  <  1  can  often  be  omitted  for  the  characteristics  of  a  motion 


over  a  relatively  long  period  of  time.  The  computation  can  begin 
with  the  boundary  layer  varying  in  a  quasi-state. 

When  a  discrete  vortex  is  newly  created,  because  the  vortex 
and  its  image  vortex  "distort”  the  local  flow  field,  the  boundary 
layer  will  be  separated  prematurely.  The  calculated  long  term 
mean  separation  angle  is  only  +67°^^  or  +  77°^^,  which  is 
smaller  than  the  widely  believed  +82°.  One  consideration  to 
minimize  this  effect  is  to  begin  counting  the  induction  effect  of 
a  recent  vortex  after  one  or  a  half  step  in  time^'^’^^.  This  is 
equivalent  to  placing  the  initial  position  of  the  recent  vortex 
at  a  distance  of  v^At  or  (1/2)vcaT  away  from  the  separation 
point,  vj  and  (1/2)vc  are  the  background  flow  field  velocity  or 
the  transfer  velocity  of  the  vortex  point.  In  conjunction  with 
viscous  vortex  nucleus  correction,  this  method  can  improve  the 
calculated  results  of  the  separation  point.  In  some  engineering 
problems  focusing  on  long  term  flow  characteristics,  in  order  to 
avoid  the  above  difficulty  in  calculating  the  separation  point, 
it  is  assumed  that  separation  occurs  at  the  place  where  the 
maximum  surface  velocity  chopped  by  This  assumption  is 

very  close  to  actual  result.  There  are  different  methods  to 
determine  the  initial  radial  position  of  a  recently  created 
vortex  by  satisfying  the  no-slip  condition  at  the  vortex  at  the 
instance  when  a  point  vortex  is  disengaged^ ^  ,  by  choosing  an 

appropriate  initial  value  from  numerical  experiments  based  on  the 
stability  of  the  vortex  motion,  and  by  using  the  initial  radius 
as  a  half  of  the  thickness  of  the  boundary  layer  to  determine  the 


initial  radial  position  of  the  vortex  when  the  viscous  vortex 

[3] 

nucleus  is  corrected  .  In  general,  in  flow  problems  involving 
the  time-dependent  separation  points,  the  initial  position  of  a 
discrete  vortex  does  not  necessarily  have  to  be  determined  by  the 
no-step  condition. 

III.  Secondary  Vortex,  Effect  of  Secondary  Vortex  and  Mechanism 
of  Vortex  Detachment 

The  numerical  modelling  of  the  discrete  vortex  method  under 
the  assumption  of  high  Reynolds  number  subcritical  conditions 
shows^-^  that  the  maximum  return  surface  flow  in  the  return 
flow  region  can  reach  the  order  of  magnitude  of  the  incoming  flow 
in  the  initial  unsteady  stage  of  the  flow  around  the  cylinder.  A 
part  of  the  rear  shear  layer  is  continuously  under  the  effect  of 
a  reverse  pressure  gradient  after  a  specific  instance. 

Theoretical  computations^-^  assumed  that  this  rear  shear  layer 
would  separate  and  create  a  secondary  vortex  in'  opposite 
direction  to  the  primary  vortex.  The  secondary  vortex  has  an 
important  effect  on  the  motion  and  evolution  of  the  primary 
vortex,  the  cancellation  of  vortex  in  the  wake,  the  distribution 
of  pressure  on  the  surface  of  the  object,  and  the  drag  and 
transverse  force  executed  on  the  cylinder^^  ’  .  Figures  1 

and  2  show  some  examples  of  the  effect  of  the  secondary  vortex. 


■1 


Figure  1  .  Effect  of  Secondary  Vortex  on  the  Evolution  of  the 
Primary  Vortex 

1.  secondary  vortex 

2.  primary  vortex 


Figure  2. 


.  ******* 
ft**** 

Effect  of  Secondary  Vortex  on  Drag  Coefficient  and 
Transverse  Force  Coefficient  C^. 

1.  results  of  computation  after  taking  secondary 
vortex  into  account 

2.  results  of  computation  by  not  considering 
secondary  vortex 
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The  vortex  motion  described  by  theoretical  calculation  * 

agrees  with  many  experimental  observations  at  high  Reynolds 
[81 

numbers  .  Therefore,  the  effect  of  rear  shear  layer  separation 
and  secondary  vortex  cannot  be  neglected.  On  the  other  hand,  in 
one  experiment  conducted  at  Reynolds  number  9500  it  was  observed  / 1 3 
that  there  was  a  pre-wake  formation  and  evolution  process 
(including  the  appearance  and  disappearance  of  a  pair  of 
secondary  vortices)  prior  to  the  formation  of  the  primary  wake. 
Afterwards,  due  to  the  instability  of  the  flow,  return  flow 
separation,  such  as  the  one  considered  in  the  theoretical 
calculation,  did  not  take  place.  As  the  calculation  of  the  rear 
shear  layer  separation  point,  the  time-dependent  variation  of  the 
separation  point  was  obtained  by  using  the  boundary  layer 
approximation  method  (Stratford  method)  under  the  assumption  of 
laminar  flow  separation  in  the  early  stage  of  the  flow  in 
reference  [7],  In  long  term  flow  situations,  there  are  only 
approximation  methods ^-51  [6]  .  Therefore,  there'  is  a  need  to 
conduct  detailed  non-contact  experiments  on  the  rear  shear  layer 
in  order  to  understand  its  special  characters  so  that  we  can 
perform  a  more  delicate  theoretical  analysis. 

There  are  different  theories  explaining  the  mechanism  vortex 
detachment  from  the  cylinder.  There  are  two  different  viewpoints 
on  the  breakdown  of  the  primary  vortex  sheet  in  its  developing 
process  due  to  the  factor  whether  the  effect  of  secondary 
vortices  should  be  considered^^  ’  .  Correspondingly, 

different  methods  were  used  to  introduce  an  asymmetric 


perturbation  to  the  numeric  calculation  in  order  to  realize 
primary  vortex  detachment.  However,  the  presence  and  development 
of  secondary  vortices  to  fracture  the  symmetric  vortex  sheet  may 
be  one  of  the  conditions  leading  the  vortex  to  begin  an 
asymmetric  motion. 

IV.  Reduction  of  Circulation  of  Vortices  in  the  Wake 
With  regard  to  the  cancellation  and  viscous  dissipation  of 
vortices  in  the  wake,  the  total  circulation  in  the  concentrated 
wake  vortex  should  be  much  less  than  the  sum  of  the  vortices 
detached  from  the  boundary  layer.  Experimentally,  it  is  only 
approximately  60%  of  the  original  value.  However,  the  calculated 
result  based  on  the  discrete  vortex  mode  1  is  around  80%  of  the 
original  value.  This  problem  has  not  been  resolved  in  the 
theoretical  model.  It  will  have  an  effect,  to  various  extents, 
on  the  separation  point,  transverse  force,  wake  length,  St 
number,  size  of  the  vortex  and  stability  region" of  vortices. 
Cancellation  of  vortices  is  caused  by  approaching  vortices  in 
opposite  directions  and  by  point  vortices  less  than  a  specific 
distance  from  the  surface  of  the  object.  There  are  three 
possible  reasons  for  opposite  direction  vortices  to  approach  one 
another  in  the  flow  behind  a  cylinder:  1)  point  vortices 
generated  by  secondary  vortices  and  boundary  layer  separation 
approach  one  another  behind  the  cylinder;  2)  distorted  or  cutoff 
primary  vortex  layer  is  carried  into  the  wake  by  the  inviscous 
flow;  and  3)  the  vortices  are  swept  apart  downstream  and  thus 
distributed  in  the  entire  wake.  Several  numerical  experiments 


were  conducted  to  study  the  extent  of  vortex  reduction  by 
different  mechanisms ^ ^  .  In  some  work,  a  correction 

factor,  less  than  1,  was  multiplied  to  the  intensity  of  recently 
created  or  secondary  vortices  in  order  to  match  the  calculated 
drag  (or  other  parameter)  with  the  experimental  results  . 

The  distribution  of  vortex  was  studied  from  various  angles.  In 
reference  [2],  the  assumption  of  point  vortex  conservation  was 
abandoned.  It  was  assumed  that  the  loss  of  each  point  vortex  in 
the  wake  is  proportional  to  its  intensity  and  its  location  in  the 
recently  separated  vortex  sheet.  Certain  principles  for 
turbulent  flow  loss  were  obtained  through  numerical  experiments. 
The  total-  circulation  was  decreased  by  approximately  207«  in  the 
calculated  result.  However,  it  is  necessary  to  further  determine 
major  causes  of  vortex  reduction  in  real  flows  and  the  patterns 
of  variation. 
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SOME  PROBLEMS  IN  DISCRETE  VORTEX  NUMERICAL 
MODELLING  ON  VORTEX  MOTION  BEHIND  A 
CIRCULAR  CYLINDER 

Ling  Guocan 

(fvritate  ol  Mechanic*,  Academia  Sinlca, Btifing) 

Abstract 

This  is  a  shortened  review  which  is  concentrated  on  four  aspects  in 
discrete  vortex  numerical  modelling  on  unsteady  flew  and  vortex  motion 
behind  a  circular  cylinder,  i.  e.  (1)  the  discrete  method  and  the  stability 
of  the  motion  of  vortex  sheet  (2)  boundary  layer  separation  and  the  nas¬ 
cent-discrete  vortices'  positions  (3)  the  secondary  vortices  problem  and 
(4)  the  redaction  of  circulation  of  vortices  in  the  wake.  Some  further 
works  are  also  suggested  in  this  note. 
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Calculation  of  Slender  Delta  Wing  with  Leading-edge  Separation 

by  Quasi-Vortex-Lattice  Method 
Xiung  Shauwen  (Beijing  Institute  of  Aeronautics  and  Astronautics) 

I.  Introduction 

The  leading-edge  separation  vortex  produced  by  a  slender 

wing  with  a  large  backward  sweep  at  a  large  attach  angle  has  an 

important  effect  on  its  aerodynamic  characteristics.  There  are 

several  theoretical  methods  to  calculate  the  aerodynamic 

properties  of  the  leading-edge  vortex.  Earlier,  there  was  the 
r  i  ] 

Smith  method  ,  which  was  based  on  a  slender  body  assumption. 
Therefore-,  the  trailing-edge  conditions  cannot  be  satisfied.  The 
leading-edge  suction  analogy  (and  its  extensions)  is  a  relatively 
simple  and  accurate  method  to  calculate  the  forces  and  moments  on 
the  entire  wing.  However,  it  does  not  provide  a  detailed  load 
distribution  .  In  reference  [3],  the  vortex  lattice  method  was 
extended  to  wings  -with  leading-edge  separation.  Several  discrete 
vortex  lines  were  used  to  represent  the  actual  leading-edge 
vortex  with  a  nucleus.  The  force  and  moment  thus  calculated  are 
more  accurate.  But,  the  load  distribution  is  not  satisfactory. 
Currently,  the  most  realistic  leading-edge  separation  model  is 
the  free  vortex  sheet  method  which  employs  a  continuous  dipole 
distribution  to  represent  the  wing  and  the  free  vortex  sheet, 
resulting  in  a  satisfactory  pressure  distribution^^.  However, 
this  method  is  more  complicated  and  requires  a  larger  computer 
and  more  computer  time. 
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Mehrota  and  Lan 


used  the  quasi-vortex-lattice  method, 


[5] 

which  is  applicable  to  flows  without  separation,  on  wings  with 

r  &  ] 

leading-edge  separation  .  This  method  divides  the  wing  along 
its  span  into  several  narrow  strips.  Each  narrow  wing  strip  is 
replaced  by  a  continuous  vortex  sheet.  The  leading-edge 
separation  vortex  is  represented  by  the  discrete  vortex  line 
added  onto  the  vortex  sheet  as  well  as  from  the  leading-edge  and 
trailing-edge.  The  special  feature  of  this  method  is  that 
boundary  conditions  can  be  met  at  the  leading-edge  through  a 
theoretical  treatment  of  singular  points  on  the  leading-edge.  In 
addition,  it  is  capable  of  taking  leading-edge  separation  into 
account.  The  latter  may  be  an  appropriate  method  to  calculate  a 
wing  strip  where  separation  of  the  inner  wing  is  different  from 
that  of  the  outer  wing.  Besides,  this  method  maintains  the 
relative  simplicity  of  the  vortex  lattice  method.  In  the 
meantime,  it  also  gives  more  satisfactory  results. 

In  this  paper,  the  theory  in  reference  [5]  is  used  to 
calculate  a  flat  delta  wing  with  complete  leading-edge  separation 
on  a  FEL/}<  C-256  computer.  Furthermore,  the  results  are  compared 
with  experimental  data. 

II.  Calculation  Method 

1.  Calculation  Model 

For  a  large  backward  sweep  thin  delta  wing,  an  attached 
vortex  sheet  represents  the  wing  surface.  The  discrete  free 
vortex  line  is  used  to  represent  the  separated  leading-edge 


36 


vortex  sheet  and  the  tail  vortex  sheet  in  order  to  simultaneously 
satisfy  the  airfoil  condition,  the  leading-edge  and  trailing-edge 
boundary  conditions,  and  the  conditions  that  no  force  is  exerted 
on  the  leading-edge  and  tail  vortices.  Figure  1  shows  the 
distribution  of  attached  vortex  elements  and  tail  vortex 
elements.  The  positions  of  attached  vortex  elements  in  the  chord 
and  span  direction  follow  the  cosine  law,  i.e. 

s.  »*i  +  -|-(l-eos  9t)  (1) 

This  paper  was  received  on  June  9,  1984. 
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where  is  the  x  coordinate  of  the  leading -edge  of  the  wing,  c 
is  the  local  chord-  length,  N  is  the  number  of  vortex  elements 
along  the  chord  direction,  b  is  the  span  length  and  M  is  the 
number  of  longitudinal  vortex  elements  parallel  to  the  plane  of 
symmetry.  The  longitudinal  vortex  elements  divide  the  half  wing 
into  (M-1)  strips  parallel  to  the  plane  of  symmetry  along  the 
span  direction.  The  vortex  density  -y  is  continuous  along  the 
chord  direction  in  every  narrow  strip.  In  the  same  narrow  strip, 
the  value  of  y  along  the  span  direction  remains  unchanged. 

The  position  of  the  control  point  (x„  ,  y  )  is  also 

cp  cp 

determined  according  to  the  law  of  cosine: 

+  i  =  0,l -AT  (3) 
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y,„=  -^l-cos-jgj,  j=l,2,-(M-l) 


where  x^  and  are  the  x-coordinate  of  the  leading-edge  and 
chord  length  at  YCpj •  When  i  =  0,  the  control  point  is  the 
leading-edge.  When  i  =  N,  the  control  point  is  in  the  trailing 
edge . 


F igure  1 . 


Distribution  of  Attached  Vortex  Elements  and  Tail 
Vortex  Elements  on  the  Wing 

1 .  control  point 
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Figure  2.  A  Typical  Leading-edge  Free  Vortex  Element 

1 .  leading-edge 

2.  attached  vortex  element 

3.  control  point 

The  initial  leading-edge  vortex  is  added  to  the  attached 
vortex  lattice.  Figure  2  shows  a  typical  leading-edge  vortex 
element  AK  on  a  narrow  strip  of  the  wing.  The  free  vortex 
element  consists  of  many  broken  vortex  segments.  The  dotted  line 
represents  the  initial  position  and  the  solid  line  indicates  the 
final  position.  A  and  K  are  a  chord  length  away  from  the 
trai ling-edge .  The  vortex  line  directions  from  infinity  to  A  and 
that  from  K  to  infinity  coincide  with  the  direction  of  the 
incoming  flow.  The  initial  position  of  each  segment  of  vortex 
line  can  be  determined  as  follows: 

(1)  Initially,  AE  is  in  the  same  direction  as  that  of 
attached  longitudinal  vortex  element.  Segment  DE  is  fixed  on 


the  wing  plane.  BC  and  CD  are  O.IC^.  In  order  to  meet  the 
trailing-edge  contact  condition,  CD  is  fixed  on  the  wing  plane 
and  BC  is  only  allowed  to  vary  in  the  x-z  plane. 

(2)  EF  is  on  the  wing  plane.  It  is  located  between  the 
leading-edge  and  the  first  attached  vortex  element.  It  is 
located  in  the  position  equivalent  to  the  first  vortex  element 
position  which  divides  the  attached  vortex  into  (N+1)  elements  in 
the  chord  direction  i.e. 

**=*,*+^-[1-cos2(7rri)] 

+  2<ATfTj] 

(3)  FH  is  on  the  wing  plane,  in  the  same  direction  as 
another  longitudinal  attached  vortex  element.  G  is  the  leading- 
edge.  In  order  to  satisfy  the  leading-edge  boundary  condition, 

GH  is  also  fixed  in  the  wing  plane. 

(4)  HK  is  rolled  up  into  the  leading-edge  vortex  on  top  of 
the  wing  by  the  end  of  the  computation.  Its  initial  position  is 
in  the  plane  parallel  to  the  plane  of  symmetry.  The  initial 
position  of  I  is: 

X 1=-XH ’  yl=yH 

z^=0 . iC^tanC 22 . 5-0 . 5a)  when  a^15° 
or  z^=0.lC^tan  a  when  a>15° 

IJ  is  in  the  flow  direction.  The  height  of  J  above  the  wing  is 
approximately  O.IC^.  JK  is  parallel  to  the  wing. 

2.  Boundary  Conditions  and  Solution 

As  described  above,  the  wing  surface  is  divided  into  (M-1) 
strips.  Each  one  has  N  unknown  attached  vortex  elements  and 


40 


one  unknown  leading-edge  vortex  element.  There  are  a  total  of 
(N+1)  (M-1)  unknowns.  There  are  (N+1)  control  points  on  each 
strip.  The  total  number  of  control  points  is  also  (N+1) (M-1). 
Boundary  conditions  are  met  at  each  control  point.  When  the 
control  point  is  not  on  the  leading-edge,  the  boundary  condition 
should  be 

)  +  C =  {  -  sin  a)  (  5a) 

If  it  is  on  the  leading-edge,  then  based  on  reference  [5],  the 
equation  is 

ZAiJM  +  iB.AiP »)  =  —  sin  a)  (5b) 


where  A. .  and  B. .  are  the  aerodynamic  effect  coefficients, 
ij  ij 

attached  vortex  element  and  leading-edge  vortex  element, 
respectively.  7  is  the  density  of  attached  vortex  elements  and  r 
is  the  condition  of  the  leading-edge  vortex  element.  a^k  is  a 
parameter  related  to  the  leading-edge  suction  coefficient.  From 


reference  [5] 


i 


where  M  is  the  Mach  number  of  the  flow,  3^  =  V'1-M^  ,  A  is  the 
leading-edge  swept  back  angle  of  the  wing,  and  Ct  is  the  cross- 
sectional  leading-edge  suction  coefficient  which  can  be  predicted 
in  advance.  If  the  leading-edge  is  completely  separated,  Ct  is 
zero  everywhere,  then  a^  =  0. 

Because  the  density  of  the  attached  vortex  elements  as  well 
as  the  intensity  and  position  of  the  leading-edge  vortex  are  not 
known, it  is  required  to  solve  the  equation  by  iteration.  The 
position  of  the  free  vortex  is  adjusted  in  the  solution  finding 
process  to  make  the  direction  of  each  vortex  segment  the  same  as 


that  of  the  flow  velocity  at  the  center.  There  are  two  standards 
, for  judging  the  convergence  of  the  solution.  One  is  the  free 
vortex  is  subject  to  the  least  amount  of  force  in  continuous 
iterations  and  the  other  is  that  the  circulation  of  the  free 
vortex  does  not  vary  by  more  than  27.  in  two  consecutive 
calculations.  If  these  two  points  are  satisfied,  then  it  is  the 
solution  we  are  looking  for. 

3.  Load  Computation  /24 

The  load  on  the  attached  vortex  segment  in  the  span 
direction  and  the  load  on  both  ends  caused  by  longitudinal 
vortices  are  first  calculated.  Then,  the  latter  is  added  to  the 
calculation  of  points  in  the  span  vortex  segments.  The  load  on  a 
point  along  the  chord  can  be  obtained  by  Fourier  series 
expansion. 


III.  Comparison  of  Computation  with  Experiments 

In  this  work  we  calculated  a  delta  wing  with  a  7 40  sweep  back,  with 
pleading  edgs^o tally  separated  and  withJ 

A  ML  =  0.  We  choose 

N=6  and  M=10.  The  length  of  the  free  vortex  segment  is  0.15C^. 

We  calculated  three  attack  angles.  In  Figures  3  to  6, 

comparisons  of  calculated  values  to  experimental  data  are  given. 

The  experimental  data  was  taken  from  reference  [7].  When  the 

attack  angle  is  10®  and  20°,  calculated  lift  agrees  with  the 

experimental  value.  When  the  attack  angle  is  30°,  the  calculated 

value  is  on  the  low  side.  On  the  curve  of  moment,  results  at  all 

these  attack  angles  agree  very  well  with  the  experimental  data. 
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As  for  the  cross-section  load  distribution  in  the  span  direction, 
calculated  results  are  generally  in  agreement  with  the 
experimental  results  near  the  trailing-edge .  Near  the  leading- 
edge,  the  calculated  peak  ACp  is  lower  than  the  experimental 
value.  Furthermore,  the  position  is  closer  to  the  root.  This  is 
due  to  the  replacement  of  the  actual  concentrated  vortex  by 
discrete  vortices.  The  aspect  ratio  of  the  wing  calculated  in 
this  work  is  relatively  small  (aspect  ratio  =  1.147).  If  the 
lattices  in  the  span  direction  can  be  increased,  the  calculation 
may  be  improved. 


Figure  3.  Comparison  of  Calculated  and  Experimental  Lift  Curve 

1.  experimental 

2.  calculated 


Figure  4.  Comparison  of  Calculated  and  Experimental  Moment 
Curve 

1 .  experimental 

2.  calculated 


(m)  X/c*»0.M#7 


Figure  5.  Distribution  of  aC^  in  the  Span  Cross-section,  a 

1.  experimental 

2.  calculated 


1.  experimental 

2.  calculated 

3.  experimental 

4.  calculated 
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A  CALCULATION  OF  SLENDER  DELTA  WING 
WITH  LEADING-EDGE  SEPARATION  BY 
QUASI- VORTEX-LATTICE  METHOD 

Xiong  Shanwen 

(Beijing  Institute  of  Aeronautics  and  Astronautics) 

Abstract 

The  Quasi- Vortex-Lattice  method  (QVLM)  which  was  used  calcula¬ 
ting  the  thin  wing  without  separation  has  been  extended  to  calculate  the 
slender  delta  wing  with  leading-edge  separation.  The  advantage  of  this 
method  is  that  the  leading-edge  boundary  condition  can  be  exactly  satis¬ 
fied.  It  can  be  used  to  predict  aerodynamic  characteristics  of  wings  having 
partial  leading-edge  separation.  A  calculation  has  been  made  here  for  a 
slender  delta  wing  with  complete  leading-edge  separation  and  the  results 
are  compared  with  the  experimental  data.  The  comparison  shows  that 
QVLM  can  give  satisfactory  or  reasonable  results. 


Calculation  of  Flow  Around  Thick  Wing  with  Separation  Vortices 
Zhu  Peiye  (Institute  of  Computer  Technology,  China  Aviation 


Institute) 

Abstract 

A  panel  method  was  developed  in  this  work  to  calculate  the 
aerodynamic  force  on  a  three-dimensional  wing  with  separation 
vortices.  The  model  used  is  simple  and  visual.  Furthermore,  it 
is  suited  for  any  arbitrary  planform  wings  with  different 
profiles.  This  method  employs  a  planar  quadrilateral  lattice  as 
well  as  another  lattice  formed  by  a  parallelogram  and  four  non- 
coplanar  triangles.  A  second  order  doublet  sheet  distribution  is 
arranged  on  each  lattice  in  order  to  ensure  a  high  accuracy. 
Examples  used  include  a  rectangular  wing  and  a  sweptback  wing. 
Calculated  results  were  compared  with  available  experimental  data 
and  other  theoretical  results.  They  were  found  in  satisfactory 
agreement . 


I.  Introduction 

A  separation  vortex  type  of  aerodynamic  distribution  is 

[  1  ] 

often  used  for  fighters  and  it  will  be  used  in  the  near  future 
In  this  layout,  the  flow  is  separated  at  the  leading  edge  of  the 
wing  tip.  A  separation  vortex  is  created.  It  flows  downstream 
from  the  upper  side  of  the  wing.  This  separation  vortex  has  a 
great  effect  on  the  pressure  distribution  on  the  wing,  resulting 


in  the  non-linear  rise  of  the  lift  coefficient  with  increasing 
angle  of  attack.  The  aerodynamic  properties  of  the  aircraft  can 
be  significantly  improved  if  this  phenomenon  can  be  effectively 
controlled  and  utilized.  The  maneuverability  will  be  improved 
dramatically.  The  typical  examples  include  the  slender  delta 
wing,  side  strip  wing  and  close  range  coupled  duck  wing  layout. 

In  recent  years,  some  work  was  done  to  calculate  the 

aerodynamic  load  on  a  wing  with  separated  vortices  and  the 

geometric  positions  of  separation  vortices.  For  instance, 

Belotserkovskii^ ^ ,  Rehback^^  and  Kandil^^,  used  a  so-called 

vortex  lattice  method,  together  with  a  relaxation  technique,  to 

calculate  the  aerodynamic  load  on  a  thin  wing  at  a  large  attack 

angle  and  the  positions  of  separation  vortices.  In  their 

calculations,  the  wing  was  replaced  by  an  arc  with  vortex  circles 

(equivalent  to  constant  doublet  distribution).  Discrete  vortex 

lines  are  used  to  represent  the  free  vortex  sheet,  in 

approximation.  The  authors  also  employed  similar  method  to 

calculate  the  aerodynamic  force  on  a  ultrathin  large  attack  wing 

r  5  Zi 

with  very  good  results  ’.  This  type  of  method  is  simple. 

But,  it  only  considers  a  thin  wing  and  cannot  calculate  the 
pressure  distribution  oh  the  wing.  Johnson  et  al  at  Boeing^ ^ 
introduced  a  panel  method  to  calculate  a  thick  three-dimensional 
wing.  In  their  model,  the  separation  vortex  consists  of  the  free 
vortex  sheet,  vortex  core  and  a  cross-section  connecting  the  two 
across  which  the  free  vortex  is  transported  to  the  vortex  core. 


A  first  order  source  distribution  was  used  for  the  middle  surface 
of  the  wing  as  well  as  for  free  vortex  sheets.  In  addition,  a 
second  order  curved  lattice  was  used  to  provide  their  method  with 
a  higher  accuracy.  A  very  good  pressure  distribution  was 
obtained.  However,  we  can  see  that  their  calculation  method  is 
more  complex  and  requires  more  computer  time. 

In  this  work,  a  method  is  developed  to  calculate  the 
aerodynamic  load  of  a  three-dimensional  thick  wing  at  large 
attack  angles.  In  this  method,  the  wing  surface  is  divided  into 
lattices.  Furthermore,  a  second  order  doublet  distribution  is 
employed.  The  free  vortex  sheet  is  represented  by  broken 
discrete  free  vortex  line  in  approximation. 


This  paper  was  received  on  March  20,  1984.  Revised  manuscript 
was  received  on  May  31. 


As  compared  to  Boeing's  method,  this  separation  vortex  model  is 
simple.  It  does  not  require  any  geometric  data  on  the  middle 
plane  of  the  wing.  There  is  also  no  need  to  calculate  the 
concentrated  source  effect  index.  It  not  only  is  applicable  to 
various  wing  profiles  but  also  has  a  satisfactorily  high 
accuracy . 


II.  Mathematical  Approach  of  the  Problem 
The  steady,  in  viscous  and  incompressible  potential  flow  of 
a  three-dimensional  thick  wing  is  considered.  It  is  assumed  that 
the  wing  tip  has  a  sharp  edge.  When  there  is  lift,  the  air  flow 
will  be  separated  at  the  wing  tip  to  create  a  separation  vortex. 
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As  we  all  know,  this  problem  satisfies  the  Laplace  equation.  The 
boundary  conditions  are  that  the  wing  surface  and  free  vortex 
sheet  Sy  are  unpenetratable ,  the  pressure  on  either  side  of  the 
free  vortex  sheet  is  equal,  the  perturbation  velocity  at  infinity 
is  zero  and  the  Kutar  condition  on  the  trailing  edge  and  line  of 
separation.  Notice  that  the  shape  of  the  free  vortex  sheet  is 
not  known  in  advance.  Therefore,  this  is  a  second  boundary  value 
problem  with  a  position  of  the  boundary  unknown.  This  unknown 
boundary  must  be  solved  by  using  an  iteration  method  as  a  part  of 
the  solution  to  the  problem.  From  Green's  formula  we  know  that 
the  solution  to  the  problem  may  be  expressed  by  a  surface 
integral  of  a  source  concentration  density  o  and  a  doublet 
density  p .  In  this  work,  only  a  doublet  surface  distribution  is 
employed.  By  taking  the  unpenetratable  boundary  condition  into 
consideration,  we  get  the  following  integral  equation 


where  indicates  that  the  Cauchy  principal  value  is  chosen  for 
the  integral,  S  =  S^+Sy ,  Vm  is  the  free  flow  velocity,  and  n  is 
the  unit  external  normal  vector.  This  is  obviously  much  simpler 
than  using  o  and  p  simultaneously.  As  for  other  methods 
involving  simultaneously  using  o  and  p  and  a  comparison  of  their 
results,  one  can  read  reference  [8]. 


III.  Numerical  Method 

1.  Approximation  of  Wing  Surface  and  Free  Vortex  Sheet 


Figure  1 .  Two  Types  of  Lattices  of  Approximate  Wing  Surface 
Area  Element 

(a)  parallelogram  lattice 

(b)  non-coplanar  plane  combination  lattice 


The  wing  surface  is  divided  into  many  small  parallelogram 
area  elements  (lattices).  Each  element  is  approximated  by  a 
parallelogram  or  by  a  non-coplanar  combination  comprised  of  four 
triangles.  Figure  la  shows  the  parallelogram  Q-j-j,  Q -j 3  *  Q-|4* 

It  is  parallel  to  two  diagonals  and  Q2Q4  °f  the  surface 

element  Q-jC^Q^Q^  on  the  wing  surface.  Furthermore,  it  is  at  the 
same  distance  away  from  them.  Figure  1b  shows  a  parallelogram 
ABCD.  The  apices  are  located  in  the  middle  of  the  lines 
connecting  the  quadrilateral  Q-|Q2^3^4*  This  parallelogram  and 
the  four  triangles,  aABQ^,  aBCC^,  ACDQ2  and  aDAQ^,  are  combined 
to  create  the  approximate  area  element  Obviously,  the 

former  is  simple  and  easy  to  calculate.  But,  the  accuracy  is  not 
very  high.  It  is  usually  used  when  the  curvature  of  the  wing 
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However,  it  can  more  accurately  approximate  the  wing  surface, 
especially  in  large  curvature  areas  such  as  the  leading  edge, 
side  and  trailing  edge.  Furthermore,  it  also  ensures  the 
continuity  of  the  geometric  shape.  In  this  work,  this  method  is 
used  to  approximate  the  wing  surface,  not  only  to  attain  certain 
accuracy  but  also  to  ensure  that  the  calculation  is  not  too 
complicated . 

According  to  Helmholtz  theory,  the  intensity  of  a  free 
vortex  sheet  remain  unchanged  downstream.  Therefore,  it  is 
possible  to  turn  a  free  vortex  sheet  into  a  discrete  vortex  lines 
of  equal  intensity.  These  vortex  lines  should  be  stream  lines. 

In  numeri-cal  computation,  they  are  approximated  by  broken  lines. 
The  last  segment  is  a  semi-infinite  line  parallel  to  the  free 
flow  velocity. 

2.  Second  Order  Doublet  Sheet  Distribution 

A  second  order  doublet  sheet  distribution  is  arranged  on 


each  lattice 


nix,y)  =  y,  +  y,x  +  y,y  y„x 1  +  y„xy  +  y„y' 


where  (x,y)  is  the  local  right  angle  coordinate  on  the  lattice. 
The  origin  is  at  the  center  of  the  lattice.  Hence,  jig  is  the 
value  of  the  doublet  distribution  at  the  center.  This  value  is 
chosen  to  be  the  discrete  unknown  of  the  problem.  Other 
coefficients  such  as  yx,  4y  •  •  •  >  ^yy  are  not  independent.  They  are 
expressed  as  a  linear  combination  of  discrete  unknowns  on  two 


neighboring  lattices 


’  1 
y>=  2  /***■=  2 


where  p^ ,  i=0 , . . . 


8  are  the  discrete  unknowns  of  the  lattice  and  8 


neighboring  lattices,  b  . ,  b  , ,...b  .  are  constants  which  can 

xi  yi  yyi 

be  ca  culated  by  the  weighted  least  square  method. 

Strictly  speaking,  the  optimal  accuracy  can  only  be  obtained 
by  matching  second  order  doublet  distribution  with  a  second  order 
curved  lattice.  But,  the  curvature  in  most  areas  on  the  wing  is 


usually  not  very  large  in  most  cases.  Because  a  non-coplanar 
combination  lattice  is  used,  in  addition  to  locally  tightening 


the  lattice,  better  accuracy  can  be  obtained.  Therefore, 
satisfactory  results  can  also  be  obtained  by  using  a  second  order 
doublet  distribution  and  a  planar  lattice. 

-  3.  Kutar  Condition 


With  regard  to  the  Kutar  condition  for  the  trailing  edge, 
there  are  many  equivalent  expressions.  In  this  method,  the 
condition  used  is  that  the  vortex  intensity  vector  on  the 
trailing  edge  is  parallel  to  the  local  velocity.  In  general,  the 
flow  velocity  at  the  trailing  edge  is  nearly  parallel  to  the  wing 
chord.  The  perpendicular  component  of  the  vortex  intensity 
vector  in  the  span  direction  is  equal  to  zero.  This  method  does 


not  require  a  new  equation  and  does  not  introduce  a  new  unknown. 


It  is  very  convenient.  “ 

4.  Formation  and  Solution  of  Linear  Algebraic  Equations 
The  surface  integral  on  the  right  hand  side  of  the  discrete 
approximation  formula  (1)  discussed  above  is  used  to  obtain  a 
series  of  linear  algebraic  equations  to  find  the  discrete  un¬ 
know  p^  by  letting  this  equation  be  satisfied  at  all  lattice 


centers  (i.e.,  control  points). 

H 

*2  o„ —  ,N 

/-»  (4) 

where  N  is  the  number  of  wing  surface  lattices  and  a„  is  the 
influencing  coefficient.  Notice  that  the  kernel  of  the  integral 
(1)  has  the  following  property 


where  S  is  an  enclosed  curve.  If  errors  in  the  discretization 
process  and  numerical  computation  can  be  neglected,  then  the 
determinant  of  the  coefficient  matrix  in  equation  (4)  |a^|  =  0- 
These  equations  have  infinite  sets  of  solutions  differing  by  a 
constant.  For  this  reason,  we  let  u^=0  and  use  the  least  square 
method  to  solve  N-1  unknowns  from  N  equations.  The  linear 
equations  are  solved  by  the  method  of  elimination.  The  method  of 
elimination  is  always  reliable  and  effective  for  equations  of 
reasonable  order. 

5.  Relaxed  Iteration  of  Free  Vortex  Line 

If  the  geometric  postion  of  the  free  vortex  is  known,  then 
the  method  described  above  can  be  used  to  find  the  doublet 
intensity  distribution  on  the  wing.  By  using  the  method  of 
superposition,  it  is  possible  to  calculate  the  velocity  at  any 
point  in  the  flow  field.  In  our  problem,  the  geometric  position 
of  the  free  vortex  is  not  known.  It  requires  a  relaxed  iteration 
technique  to  simultaneously  determine  the  doublet  intensity  and 
the  geometric  position  of  the  free  vortex.  This  relaxed  itera¬ 
tion  process  can  be  simply  described  in  the  following. 


The  doublet  intensity  distribution  is  calculated  by  using 
the  above  method  with  a  given  initial  free  vortex  position,  i.e., 
given  free  vortex  lines  from  the  trailing  edge  of  the  wing  as 
well  as  from  the  separation  line.  In  the  flow  field  defined  by 
this  doublet  intensity  distribution,  the  position  of  the  free 
vortex  line  is  corrected  based  on  the  condition  that  a  free 
vortex  line  is  a  stream  line  in  order  to  make  it  closer  to  the 
stream  line.  The  calculation  is  repeated  by  using  the  corrected 
free  vortex  line  as  the  initial  value  until  the  maximum 
correction  of  two  consecutive  iterations  is  less  than  a  specific 
accuracy  requirement. 


IV.  Examples 

1.  Rectangular  Wing 

A  symmetric  rectangular  wing  with  an  aspect  ratio  2,  the 
Boeing  TR17,  is  under  consideration.  The  distribution  of  wing 
thickness  along  the  span  direction  is 

d\y)  =i.  11  cos-7-  ,Jt  ^£[0,11  (5) 

The  half  wing  span  is  divided  into  4  lattices  along  the  span  and 
12  lattices  in  the  chord  direction  according  to  the  law  of  cosine 
to  allow  the  leading  edge  to  have  a  tighter  mesh.  Thus,  the  half 
wing  span  has  96  lattices.  Figure  2  shows  the  calculated 
separation  vortex. 
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Figure  2. 


Side  Separation  Vortex  of  a  Symmetric  Rectangular 
Wing 

(a)  separation  vortex  in  space  o  =  19° 

(b)  cross-section  of  separation  vortex 
across  the  x  =  const,  plane,  a  ~  10° 


Figure  3.  Lift  Coefficient  and  Pitch  Moment  Coefficient  of  a 
Rectangular  Wing 

1 .  This  work 

2.  lift  line  theory 

3.  vortex  lattice  method 

4.  this  method  without  separation  vortex 
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Figure  2a  shows  the  separation  vortex  at  an  attack  angle  a  = 
Figure  2b  shows  the  cross-section  of  the  separation  vortex  across 
a  plane  x  =  const,  in  the  span  direction  when  a  =  10°.  From 
these  two  figures  one  can  see  the  curling  of  the  separation 
vortex  and  the  spatial  position  of  the  free  vortex  line. 

Figure  3  is  a  comparison  of  total  lift  coefficient  and  total 

pitch  moment  coefficient  calculated  by  using  this  method  to  those 

obtained  with  other  methods.  In  the  figure,  curve  1  is  obtained 

using  the  method,  curve  2  is  based  on  the  lift  line  theory,  curve 

[9] 

.3  is  by  using  vortex  lattice  method  ,  and  curve  4  is  also  based 
on  this  method  but  without  considering  separation  vortices.  When 
the  attack  angle  is  small,  all  methods  are  comparable.  With 
increasing  attack  angle,  this  method  results  in  larger  non-linear 
lift  and  pitch  moment. 

In  order  to  demonstrate  the  adaptability  of  this  method  to 
various  wing  profiles,  it  was  used  to  calculate  the  same  wing 
with  an  asymmetric  profile,  the  E-603  wing.  The 'results  are 
satisfactory.  Figure  4a  is  the  separation  vortex  diagram 
calculated  at  o  *  22°.  An  interesting  result  was  obtained  when 
we  calculated  the  situation  that  a  =  10°.  In  this  example,  the 

wing  tip  is  not  straight.  Instead,  it  is  the  middle  arc  of  the  E-603 
profile.  At  the  leading  edge,  this  middle  arc  has  a  local  angle  of 
attack  p=20°  and  since  c£=10o,  a  local  negative  angle  of  attack  of 
approximately  -10°  occurs  at  the  leading  edge.  Therefore,  the  separation 
begins  curling  outward  and  does  not  curl  inward  as  under  normal  conditions* 

T'  e  cross-sectional  view  of  the  separation  vortex  on  the  x=const. 


plane  shown  in  Figure  4b  clearly  illustrates  this  point.  Each 
cross-sectional  curve  begins  with  an  outward  curl  and  then  curls 
inward  (as  compared  to  Figure  2b).  This  shows  that  our  method 
can  be  used  to  calculate  complicated  separation  vortices. 


Figure  4.  Side-fringe  Separation  Vortex  of  a  Rectangular  Wing 
With  Asymmetric  Profile 

(a)  spatial  diagram  of  separation  vortex,  a.  =  22° 

(b)  cross-section  of  separation  vortex  at  x  =  const., 
a  =  10° 


Figure  5.  Wing  Top  Separation  for  a  45°  Sweepback  Wing  With 
Aspect  Ratio  2  at  a  =  10° 

(a)  spatial  diagram  of  separation  vortex 

(b)  cross-section  of  separation  vortex  on  x  =  const, 
plane 


2.  Sweepback-  Wing 

Let  us  consider  an  equal  chord  length,  45°  sweepback  wing 
with  an  aspect  ratio  of  2.  This  wing  has  the  same  profile  and 
span  thickness  distribution  as  those  for  the  above  rectangular 
wing.  Figure  5a  shows  the  spatial  shape  of  the  calculated 
separation  vortex  at  a  =  10°.  Figure  5b  shows  the  cross-section 
of  this  vortex  on  the  x  =  const,  plane.  We  can  see  that  this 
method  is  also  capable  of  calculating  the  geometric  shape  of  the 
wing  top  separation  vortex  of  a  sweepback  wing. 


Figure  6.  Total  Lift  Coefficient  of  the  Sweepback  Wing 

1.  this  method 

2.  linearization  theory 
o  experimental 

Figure  6  shows  the  calculated  total  lift  coefficient  of  the 
sweepback  wing  (curve  1)  and  compared  it  with  results  obtained 
using  other  methods,  as  well  as  with  experimental  data.  Curve  2 
is  based  on  a  linearization  theory,  curve  3  is  based  on  Gers ten's 
theory,  and  curve  4  is  based  on  Bollay's  theory.  These 
theoretical  results  and  "the  experimental  data  are  directly 
brought  over  from  Reference  [10].  At  small  attack  angles,  all 
theoretical  results  are  in  agreement  with  the  experimental  data. 
At  medium  attack  angles,  the  results  obtained  with  this  method 
fit  the  data  better.  As  the  attack  angle  further  increases,  such 
as  greater  than  20°,  all  theoretical  results  deviate 


significantly  from  the  experimental  data.  This  may  be  due  to  the 
breakage  of  the  separation  vortex  over  the  wing  in  the  experiment 
which  leads  to  a  rapid  drop  of  the  lift. 

V.  Conclusion 

In  this  work,  a  panel  method  was  developed  to  calculate  the 
flow  around  a  three-dimensional  thick  wing  with  separation 
vortices.  This  method  employs  a  simple  model  and  is  applicable 
to  any  arbitrary  planar  wing  of  various  profiles  as  long  as  the 
wing  has  sharp  fringes.  Examples  shown  in  this  paper  indicate 
that  the  accuracy  of  calculated  results  is  satisfactory. 

-  The  method  not  only  can  be  used  to  calculate  the  flow  around 
a  thick  wing  with  separation  vortices,  but  also  can  be  extended 
to  complicated  situations  such  as  leading  edge  separation  vortex 
and  leading  and  trailing  edge  interference  of  a  slender  wing. 

For  a  smooth  body  with  symmetric  separation  vortices,  such  as  the 
fuselage  and  missile  head,  as  long  as  the  position  of  the 
separation  line  is  known,  this  method  can  also  be  used  to  perform 
the  calculation.  By  using  the  Planck-Glow  theory,  this  method 
can  easily  be  extended  to  calculate  a  subsonic  subcritical  flow. 

Appendix 

This  work  was  completed  by  the  author  at  the  Institute  of 
Mechanics  at  Stogadtt  University  in  West  Germany.  All  the 
computations  were  carried  out  in  the  university  computer  center. 
The  author  wishes  to  express  his  gratitude  to  the  director  of  the 
institute  and  his  advisor  Professor  R.  Eppler  for  their  assistance. 
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Zhu  Peiye 

(Institute  of  Computing  Technology,  CAE) 

Abstract 

The  present  paper  has  developed  a  panel  method  predicting  the  non¬ 
linear  aerodynamic  loads  on  thick  wings  with  separation  vortices.  The 
model  used  is  simple  and  visual.  The  method  can  be  used  for  arbitrary 
planform  wings  with  different  profiles.  Planar  quadrilateral  panel  and 
panel  that  consists  of  a  parallelogram  and  four  triangles  are  used.  In 
order  to  obtain  a  high  accuracy,  the  wing  is  represented  by  piecewise 
continuous  quadratic  doublet  sheet  distributions.  The  aerodynamic  loads 
on  rectangular  wing  and  sweepback  wing  are  computed.  They  agree  well 
with  experimental  tests  and  other  theories. 
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The  Split  Coefficient  Matrix  Difference  Method  for  Supersonic 

Three-dimensional  Flow 


/  41 


Zhang  Lumin  and  Shan  Xiaonan 
(China  Aerodynamic  Research  and  Development  Center) 

Abstract 

In  this  paper,  a  split  coefficient  matrix  (SCM)  difference 
method  is  introduced  to  calculate  the  inviscid  steady  flow  over  a 
non-symmetr ic  three-dimensional  body.  This  method  is  based  on 
the  mathematical  theory  of  characteristics.  In  the  split 
coefficient  difference  method,  the  coefficients  are  split 
according  to  the  slope  of  the  characteristic  line.  These 
coefficients  are  multiplied  by  an  appropriate  unidirectional 
difference.  The  forward  difference  is  related  to  the  negative 
characteristic  value.  The  backward  difference  is  related  to  the 
positive  characteristic  value. 

A  blunt  spherical  cone  is  used  as  an  example  in  this  paper. 
Furthermore,  we  compared  this  method  with  reference  [2]  and  found 
that  the  SCM  can  still  maintain  a  high  accuracy  with  few  meshes. 

I.  Basic  Equations 

In  the  cylindrical  coordinate  0-zr<p,  let  us  assume  that  the 
body  surface  and  shock  wave  surface  are  respectively  expressed 

3S  '  r*  =  rtU,?>)  r,  =  r,(z,<p) 
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Let  us  introduce  the  following  transformation 

2  =  2, 

where  r,  z,<p  represent  the  radial,  axial  and  circumferential 
coordinate  of  the  cylindrical  coordinate,  respectively,  r^  and  rg 
are  the  body  and  shock  wave  radial  coordinates. 

By  making  the  variables  dimensionless,  i.e.,  dividing  p  and  p 
by  flow  parameters  pB  and  pm,  dividing  velocity  by  /p^/p^,  and  by 
dividing  lengths  such  as  z  and  r  by  the  radius  of  the  blunt  RN , 
the  three-dimensional  inviscid  steady  Euler  equations  are 

Ad.  +  Bd  |  +C</#+Z>**0 

i.c,  (1) 

C.-yMl+JL. 

where  u,  v,  w  are  the  axial,  radial  and  circumferential  velocity 
components,  7  is  the  specific  heat,  p  is  pressure  and  p  is 
density. 

d=(u,v,w,p)^  pa*  0  0  «• 

P®  0.0  1 

D  0  0  . 

■  #...»•  P*  O' 


This  paper  was  received  on  January  9,  1984.  Revised  manuscript 
was  received  on  May  22. 


*?■ 

i'll 

i 


n 


pa’  l. 

pa’  £ 

PU 

0 

0 

P# 

o  i. 

o  I. 

PU 

-b"* 


I_1_PU) 

*■  ■ 


—pw 


—pa’  v  ! 

r 


D=  1  ,  i 

pw1  i 


—  ,o:-a’  1 


t7=«|,  +  y|,  +— f* 

‘"7^ 


rbz’  rb<p’  rsz  anc*  rsip  are  ^*rst  derivatives  of  or  rg  with 
respect  to  the  subscript.  E  ,  £  and  5  are  first  derivatives  of 

L  Z  <p 

with  respect  to  the  subscript. 

II.  Split  Coefficient  Equation 

The  Euler  equation  (1)  is  hypobolic.  Hence,  there  exists  a 

similarity  transformation  which  can  be  written  as 

d.  +  A-'  TA,*T"  Adt+A-'SAcjS-'Ad'  +  A- 
B-A-'B-A-'TAmmT-'A 
Let  C  -A-'C-A-'SAtjS-'A 


warns 
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Based  on  the  sign  of  the  characteristic  value,  the  diagonal 
matrix  A  can  be  divided  into  two  parts;  A+=  A+  A  where  A+  is 
only  x>0  and  A_  for  x<0. 

Therefore,  the  split  coefficient  matrix  can  be  separated  as 

§  -  A"TAba.  T-  A+A *»  TAu-  T~*  A-j}  ++B  - 
C-A-'  SAca.+S‘,A+A-'SACa-S-'  A-c.+C- 
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The  split  coefficient  equation  for  the  three-dimensional  steady 


Euler  equation  is: 

*■+ $ *dt,+ B  € .dif+A-'D-Q  '  ( 

The  subscript  b  represents  the  backward  difference  and  f 
represents  the  forward  difference.  The  above  equation  is  an 
interior  point  calculation  equation. 


III.  Boundary  Treatment 
(1)  Body  Surface  Boundary  (5  =  0) 

Based  on  the  steady  Euler  equation  (2),  the  characteristic  / 44 
compatibility  relation  is  established  as  follows: 

T-'  dt  +  T-lBdi,  +  T  '(C  +  c  d;)  +  T  'A  'D=  0 

On  the  body  surface,  the  impermeable  condition  should  be 
satisfied : 

(7  =  0  A}.,*  0 

B 

The  negative  characteristic  value  X^  should  be  taken  according  to 
the  physical  significance.  Using  9p/az  as  an  example: 

The  split  coefficient  expansion  is: 

P.=  ~  Ui  Pt-i-g.)-  ~j’-  |^r»,  itf  _  + 

A  +  nr,„  4  w  (  r>„  -  r±’-L‘i- 

.  r*  7  L  \  rt 

where 


|ts.  rt. 


9  ~  (9u  9i>  9:t  9»)T 

9-<C*d»  +  Cd9,)  +  A-'D 
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The  subscripts  b  and  f  represent  the  backward  and  forward 
difference,  respectively. 


(2)  Shock  Wave  Boundary  (£  =  1) 

On  the  shock  wave  surface,  the  Rankine-Hugoniot  relation  is 


used : 


2  2 


P  =  p ) 

r<i=  .  -L  -  r,./r,)  +  w..x 


xv'd'.-tu. r„/r,y  +  (.ul-v.*)  (l  +  rj./rj  j } 

■'-».  +  -..(l  +  <■;.+ r:./r‘ 

w  *■  w. — (v — v.)r,  ,/r, 

v  is  the  normal  velocity  in  front  of  the  wave,  u  ,  v  and  w 

n®  00  7  00  00 

are  the  incoming  flow  velocity  components. 


.y//: 

ry  y  Mm^/l  —  c®»*  a  cm'fi-cos  Cf> -o) 

■’  _  j-.  -  >  .  .1 

«-t*-*(ti0/si*a) 


a  and  B  are  the  attach  angle  and  slip  angle,  respectively. 

As  for  pressure  behind  the  shock  wave  p,  the  pressure 
characteristic  compatibility  relation  is  established  by  using 
equation  (3).  In  this  case,  we  should  take  the  positive 
characteristic  value  and  d^  is  changed  to  d^* 

IV.  Difference  Scheme 

The  MacCormack  second  order  explicit  scheme  is  used  in  this 
r  1  ] 

work  .  In  order  to  adapt  to  the  special  feature  of  the  split 


coefficient  method  and  to  ensure  the  second  order  accuracy,  a  two 

point  and  a  three  point  formula  are  used. 

A >P  - 

±<P  A<P  . 

V>j>  _  P,~P,-x 
A<P  A  <t> 

3,  p  _ 

A  <P  A<£ 

\p  _  2p,  —  3p,-i  +  pi-i 
\<t>  A<p 

If  Af  and  6^  are  used  in  predicting  step,  then  and  6^  should 
be  used  in  the  correcting  step. 


V.  Stability  Condition 

The  step  Az  must  be  limited  by  the  stability  condition.  For 
a  MacCormack  explicit  difference  equation,  stability  is  a 
necessary  condition.  With  respect  to  any  point  in  the  area  where 
solution  is  sought,  the  development  region  of  the  differential 
equation  must  be  totally  included  in  the  dependent  area  of  the 
difference  equation,  i.e. 


where 


a  = 


u  —  O' 


-M-  ->  max  {(C7,),.i,  Ccr, > , . » } 

Az 

i  <  t  « 1 1 


•roV  (p£,  +  u>£,/r)'+ (m'-o’)(£J -f  £{/r')  } 


r(u‘-a‘)  A  0 


+  Ov^  u’  +  w’-a’  } 


{ ! *  [  u*‘+vt’+ -r(£»+il)] -0’ | + 

VH- $)]i+  ]  1 
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Therefore 


/  4 


A *=  n&£/maz{(o1),.t,  (a2),.„ 

i<i </ i 


VI.  Calculated  Results  and  Analysis 

In  this  paper,  the  equation  to  calculate  the  supersonic  flow 

of  a  non-symmetric  (r^  ^  0)  re-entry  warhead  with  a.  and  &  is 

given.  In  order  to  explain  the  special  feature  of  this  method  we 

calculated  the  flow  fields  and  aerodynamic  forces  of  a  9° 

spherical  cone  at  M  =  5,  7,  9  and  a  =  0° ,  3°,  6°,  9°,  10°.  The 

results  are  compared  to  those  obtained  with  the  conventional 

ro] 

difference  method  (see  Tables  1  and  2). 


Table  1.  Comparison  of  Pressure  Coefficients  C  at  M  *  5, 
a  =  0°  p 

zCaxial  position)  difference(27x13)  dif ferenceC 1 1x7)  SCM(11x7) 

30  0.0575  0.0548  0.0586 


36 


0.0588 


0.0563 


0.0598 


Table  2.  Comparison  of  Pressure  Center  (Cm/CM)  at  M  =  5, 

30  rn  ™ 

zCaxial  position)  difference(27x13)  dif f erenceC 1 1x7)  SCM(11x7) 

30  0.7165  0.7200  0.7131 

36  0.7114  0.7145  0.7086 

From  the  above  tables  we  know  that  the  SCM  method  can  still 
maintain  a  high  accuracy  with  few  meshes.  Due  to  the  fact  that 
the  split  coefficient  method  has  good  physical  significance,  it 
is  more  appropriate  than  the  conventional  difference  method  when 
the  transverse  flow  gradient  is  large. 
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THE  SPLIT-COEFFICIENT  MATRIX  METHOD  FOR 
SUPERSONIC  THREE  DIMENSIONAL  PLOW 

"  7  “ 

*  r  * 

Zhang  Lumin  .  Shan  Xiaonan 
(C kina  Aerodynamic  Research  and  Development  Centre > 

Abstract 

The  Split  Coefficient  Matrix  (SCM)  finite  difference  method  for 
solving  inviscid  steady  supersonic  flow  over  non-symmetrical  body  is 
presented.  This  method  is  based  on  the  mathematical  theory  of  characte¬ 
ristics.  In  the  SCM  approach,  these  coefficients  are  split  according  to  the 
sign  of  the  characteristic  slopes.  The  split  coefficients  are  multiplied  by 
appropriate  one-sided.  Forward  differences  are  associated  with  negative 
characteristic  slopes.  While  backward  differences  are  associated  with 
positive  slope  values. 

The  numerical  example  of  the  blunt  sphere-cones  have  been  worked 
in  this  paper  and  compared  with  ref.  2  to  demonstrate  good  accuracy  of 
SCM  in  rare  mesh. 


Investigation  on  Some  Characteristics  of  the  FLIC  Method  and  Its  /48 
Application  to  Calculation  of  Pitot  Pressure  of  Dusty-gas  Shock 

Tube  Flow 
Du  Xixin 

(China  Aerodynamic  Research  and  Development  Center) 

Abstract 

In  this  paper,  a  trial  method  is  used  to  explore  some  of  the 
features  of  the  FLIC  method.  The  original  mass  flow  calculation 
method  is  modified.  This  modified  FLIC  method  is  used  to 
calculate  the  Pitot  pressure  of  a  dusty-gas  flow.  The  result 
reveals  the  variation  of  pitot  pressure  with  loading  ratio  and 
time.  The  stationary  part  of  the  Pitot  pressure  agrees  with  the 
analytical  solution  of  the  effective  flow  in  the  equilibrium 
region. 


I.  Introduction 

The  FLIC  method  have  many  attractive  features.  (1)  As 

compared  to  the  point  mass  lattice  method,  it  does  not  require 

many  point  masses.  For  a  given  problem,  this  will  significantly 

reduce  memory  storage  and  computer  time.  Furthermore,  the 

fluctuation  which  is  characteristic  to  the  point  mass  lattice 

method  can  be  avoided.  (2)  It  can  be  successfully  applied  to  the 

ri] 

calculation  of  two-dimensional  flow  problems  .  However,  just 
as  other  difference  methods,  an  artificial  viscosity  term  must  be 
introduced  to  the  FLIC  method  to  improve  the  stability  of  the 
calculation.  Just  as  reference  [23  pointed  out  that  the 


artificial  viscosity  term  will  blur  the  shock  wave  and  contact 
surface.  For  a  dusty-gas  flow,  details  of  the  flow 
characteristics  caused  by  the  presence  of  dusty  particles  will  be 
smeared.  The  flow  pattern  is  distorted.  Even,  the  reality  may 
be  destroyed.  To  this  end,  an  attempt  is  made  in  this  work  to 
minimize  the  appearance  of  discontinuous  surfaces  in  FLIC 
calculations . 

As  we  know  very  well  that  dusty-gas  shock  tube  flow  problems 

have  very  important  practical  significance.  Until  now,  this 

[3-7] 

problem  has  already  been  widely  studied  .  Therefore,  some 

basic  phenomena  in  a  dusty-gas  shock  tube  flow,  such  as  the 
effect  of  particles  on  the  structure  of  the  shock  wave  and 
contact  surface,  the  presence  of  a  non-equilibrium  region  behind 
the  shock  wave,  and  the  formation  of  an  effective  flow  in  the 
equilibrium  region  to  follow,  are  thoroughly  understood. 

However,  studies  on  the  characteristics  and  variation  of  the 
effective  flow  in  the  equilibrium  region  and  two-dimensional 
dusty-gas  flow  problem  are  still  inadequate.  Therefore,  there  is 
a  need  to  further  investigate  these  problems. 

This  paper  was  received- on  January  11,  1984.  Revised  manuscript 
was  received  on  June  14. 

II.  The  FLIC  Method 
1.  Brief  Introduction 

A  detailed  introduction  of  the  FLIC  method  can  be  formed  in 


reference  [1].  Here,  we  are  only  going  to  explain  its  major 
characteristics.  In  order  to  find  the  solution,  a  volume 


including  the  fluid  is  divided  into  several  lattices,  as  shown  in 
Figure  1.  The  center  coordinate  of  a  typical  lattice  is  i,  j, 
which  is  determined  by  the  position  (i-1/2)6x  and  (j-1/2)6y.  A 
finite  difference  approximation  method  for  the  equation  of  motion 
of  the  fluid  is  used  to  bring  the  quantity  in  each  lattice 
forward  with  time.  Each  computation  cycle  consists  of  two  steps. 
First,  the  median  values  of  velocity  and  specific  internal  energy 
are  calculated  by  considering  the  accelerating  effect  of  the 
pressure  gradient  alone.  This  means  that  there  is  no  mass  flow 
across  the  boundary  of  the  lattice.  The  second  step  includes  the 
calculation  of  the  migration  effect.  It  is  assumed  that  the  mass 
flow  across  the  lattice  boundary  is  directly  proportional  to  the 
density  of  flow  out  of  the  lattice.  This  is  the  so-called  "flow 
out  of  lattice"  difference  method.  The  advantage  of  this  method 
is  that  the  calculation  is  very  stable  in  the  low  velocity 
region.  Furthermore,  the  possibility  that  a  lattice  is  empty  is 
eliminated.  Therefore,  the  final  values  of  density,  velocity  and 
specific  internal  energy  can  be  obtained  from  the  conservation  of 
mass,  momentum  and  energy.  As  time  progresses,  these 
calculations  are  repeated  for  every  time  increment. 


Figure  1.  Schematic  Diagram  of  the  Calculation  Lattice 


2.  Exploration  of  Some  Features 

Two  different  schemes  are  used  to  calculate  an  ideal  shock 
tube  flow.  Figure  2  shows  the  result  obtained  using  the  original 
method  in  reference  [1],  where  the  coefficient  for  the  artificial 
viscosity  term  B  =  0.  In  the  original  method,  the  equation  to 
calculate  the  mass  flow  across  the  lattice  boundary  is  as 
follows:  if  u1?  i  .  >  0 

a*:  (D 


If  “i  +%:  j ’  <  °»  then  - 


$ *+i p' * 


(2) 
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Figure  2.  Solution  to  an  Ideal  Shock  Tube  Wave  (p^/p^=10.  B=0) 

•  a.  density  distribution 
b.  velocity  distribution 

where  aM  is  the  total  mass  flow  across  the  right  lattice  boundary 
within  a  time  increment  6t,  is  the  area  of  the  right 

lattice  boundary,  P^,j  or  1 » j  so-called  density  of  the 

flow  into  our  out  of  the  lattice,  u  is  the  x-component  of  the 
median  velocity,  and  the  superscript  n  represents  the  cycle 
number  of  the  calculation,  u1?  ,  .  =  (1/2)  (u?  .  +uV  -  ,  ). 

From  Figure^,  we  can  see  that  when  B  =  0  there  are  obvious 
overshoots  and  fluctuations  for  parameters  such  as  velocity  and 
density  behind  the  shock  wave.  Furthermore,  the  contact  surface 
is  blurred.  In  the  FLIC  method,  the  calculation  of  the  mass  flow 
is  extremely  important.  In  order  to  overcome  these  difficulties, 
the  original  method  to  calculate  the  mass  flow  is  modified. 

Figure  3  shows  the  result  of  the  modified  mass  flow  calculation 


method . 


If  u*  . 
i+h;  j 


The  mass 
>0 


flow  calculation  equation 


is : 


<0,  th6n 


AA/*.$;,  =  it 

AAf!„}.,=  p’.i.,  u 


(3) 

(4) 


'  From  Figure  3  one  can  see  that  the  shock  wave  front  and  the 
contact  surface  are  very  clear.  Furthermore,  there  is  no 
overshoot  and  oscillation  behind  the  shock  wave  front.  However, 
with  increasing  calculation  time,  it  becomes  unstable,  as  shown 
in  Figure  4. 


a.  density  distribution 

b.  velocity  distribution 

3.  Overcoming  the  Difficulties 

In  order  to  overcome  this  difficulty,  these  two  methods  are 


used  alternatively  to  stabilize  the  calculation.  In  addition, 
the  overshoot  and  oscillation  become  very  small.  The 


discontinuous  front,  however,  remains  very  clear,  as  shown  in 
Figure  5. 


Figure  4.  Solution  to  an  Ideal  Shock  Tube  (p< /p1  =  10, 
'  B  =  0) 

1.  density  distribution 


Figure  5.  Solution  to  an  Ideal  Shock  Tube 
B  =  0) 

1.  density  distribution 
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III.  Calculation  of  Pitot  Pressure  of  Unsteady  Dusty-gas  Shock  / 5 1 


Tube  Flow 

Pitot  pressure  is  a  very  important  aerodynamic  parameter. 
Methods  to  calculate  and  measure  it  have  very  important  practical 
and  theoretical  significance.  For  dusty-gas  flows,  the  study  on 
Pitot  pressure  is  not  enough.  To  this  end,  the  aforementioned 
modified  FLIC  method  is  used  to  calculate  the  pitot  pressure  of 
an  unsteady  dusty-gas  flow. 

1.  Basic  Equations 


Let  us  assume  that  a  dusty-gas  mixture  flow  through  a  tube. 
In  order  to  determine  the  equation  of  motion  for  this  mixture,  it 
is  necessary  to  make  some  assumptions.  The  gas  is  a  total  gas. 
With  the  exception  of  interaction  with  solid  particles,  its 
viscosity  and  thermal  conductivity  can  be  neglected.  The 
particles  are  uniform  spheres  which  are  incompressible.  There  is 
no  variation  of  mass.  Their  number  must  be  large  enough  so  that 
they  can  be  treated  continuously.  However,  there  must  be  so  few 
that  that  their  total  volume  can  be  neglected.  Other  detailed 
assumptions  can  be  found  in  reference  [6]. 

Let  p,  p,  T,  u  and  v  correspond  to  the  pressure,  density, 

S  8 

temperature  and  x  and  ycomponents  of  the  velocity  of  the  gas, 
repsectively .  Let  o,  (£D,  u^  and  vp  correspond  to  the  mass 
concentration,  temperature,  and  x  and  y  components  of  the 
velocity  of  the  particle,  respectively.  Under  the  above 
assumptions,  the  equations  of  mass,  momentum  and  energy 

c  ,  .  u  .  fff) 

conservation  of  the  gas  and  the  dust  are: 


Continuity  equations: 


Momentum  equations: 
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Energy  equations : 
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where  U  =  (u^  +  v^)^  I  =  C.T,  I  =  C  (3,  m  is  the  mass  of  a 
g  g  g  g  V  9  m 

single  particle,  and  and  represent  the  x  and  y  components 
of  the  drag  D  on  the  sphere  by  the  flow.  Q  represents  the  heat  / 52 
exchange  between  the  particle  and  the  gas.  The  equation  of  state 
of  the  gas  can  be  expressed  as: 


p  =  (-y-DlgP 


(13) 


From  these  equations,  the  mutual  interacting  factors  between 


the  particle  and  the  gas  are  the  drag  D  and  the  heat  exchange 


rate  Q.  In  order  to  obtain  a  set  of  closed  equations,  the 
dependence  of  D  and  Q  on  the  flow  parameters  must  be  given.  In 
this  work,  the  following  drag  equation  is  used: 


D  =  J^d-  p(U,-Ut)\U,-U,\CD  (14) 

8 

where 

Cc=0. 48-28  Re  '  "  (15) 

The  heat  exchange  formula  is: 


Q=xd  (iC,  Pr  (T-O)Nu  (16) 


where 


Vm  =  2.0  +  0.6  Pr Re"' 


(17) 


where  Cp  is  the  isobaric  specific  heat  of  the  gas,  p.  is  the 
viscosity  coefficient  of  the  gas,  and  Re  is  the  Reynolds  number 
calculated  based  on  the  diameter  d  of  the  particle  and  the 
relative  velocity  of  the  particle: 

Re=p,\U,-U,\d/n  (-jg) 


Pr  is  the  Prandt  number 

Pr=nC,!K  (-19) 

where  K  is  the  thermal  conductivity  of  the  gas. 

Because  the  surrounding  environment  of  the  particle  may 
experience  very  fast  temperature  changes,  the  temperature 
dependence  of  the  viscosity  coefficient  /j  and  thermal 
conductivity  K  must  be  considered.  The  viscosity  coefficient  of 
air  can  be  expressed  as: 


where  pr  is  the  viscosity  coefficient  at  the  reference 

temperature  T  .  Because  the  temperature  dependence  of  K  is 

similar  to  that  of  u  and  Cp  is  a  constant,  the  Prandt  constant  is 

still  a  constant  Pr  =  0.75.  In  this  calculation,  the  particle  is 

3 

a  small  glass  ball.  The  density  of  the  material  is  2.5g/cm  and 
the  diameter  is  10  micron. 

2.  Effect  Flow  and  Its  Pitot  Pressure  in  the  Equilibrium 
Region 

From  references  C 6 3  and  [7]  we  know  that  the  particles  in 
the  flow  reach  the  limiting  velocity  and  temperature  of  the  gas 
after  a  relaxation  region  behind  the  shock  wave,  i.e.,  the  dust- 
gas  -mixture  reaches  an  equilibrium.  The  mixture,  in  equilibrium, 
may  be  considered  as  a  new  total  gas^*^  -  an  effective  gas.  The 
effective  parameters  such  as  specific  heat  and  soundspeed  of  this 
effective  gas  had  been  derived  and  explained  in  reference  [6]. 

In  order  to  facilitate  the  following  description,  we  are 
introducing  two  important  parameters  here:  (I)1  frozen  Mach 
number  M^  which  is  defined  as  the  ratio  of  the  flow  to  the  local 
gas  sonic  speed  and  (2)  equilibrium  Mach  number  Mg  which  is 
defined  the  ratio  of  the  flow  to  the  local  effective  sonic  speed 
of  the  mixture. 

For  this  type  of  effective  flow  in  equilibrium,  the  Pitot 
pressure  can  be  expressed  as:  When  M  <  1 


/ 


When  M  >  1 
e 


IV.  Results  and  Discussion 

The  modified  FLIC  method  is  used  to  solve  the  partial 
differential  equations  (5)  -  (12)  for  a  dusty-gas  shock  tube 
flow.  Furthermore,  the  variation  of  Pitot  pressure  with  time  is 
calculated  at  various  loading  ratios. 

1.  Dust-gas  Shock  Tube  Flow 

Figure  6  shows  the  velocity  distribution  (using y  P^/p*  to 
make  it  dimensionless)  and  density  distribution  (using  po  to  make 
it  dimensionless)  of  a  dusty-gas  shock  tube  flow.  The  distance  x 
is  made  dimensionless  by  using  1  =  (4/3)  (pp/poo)d.  p^  is  the 
density  of  the  material.  The  high  pressure  section  of  the  shock 
tube  is  filled  with  pure  air  and  low  pressure  section  is  filled 
with  the  dust-gas  mixture.  The  diaphragm  pressure  ratio  P^/P-j  = 
10.  From  Figure  6  we  can  see  that  the  presence  of  particles 
makes  the  shock  wave  decay.  It  is  followed  by  a  non-equilibrium 
region  where  the  particles  exchange  momentum  and  energy  with  the 
gas.  When  the  velocity  and  temperature  of  the  particle  reach  the 
corresponding  values  of  the  gas,  the  flow  is  in  equilibrium. 

From  the  density  distribution  we  can  see  that  the  decaying 
contact  surface  has  a  fixed  area.  The  calculation  grasped  the 
major  features  of  the  dusty-gas  shock  tube  flow.  In  this  sense, 
the  modified  FLIC  method  is  reliable. 
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2.  Pitot  Pressure 


Figure  7  shows  the  variation  of  Pitot  pressure  with  time. 

The  loading  ratio  a  =  1.0,  M  =3.6  and  d  =10  micron.  One  can 

s  p 

see  that  a  pressure  jump  is  created  when  the  shock  wave  sweeps 
across  the  pressure  probe.  After  an  unsteady  process,  Pitot 
pressure  maintains  at  a  constant  level.  The  flow  between  the 
shock  wave  front  and  the  end  of  the  probe  is  a  non-equilibrium 
flow.  Therefore,  the  Pitot  pressure  measured  is  the  non¬ 
equilibrium  Pitot  pressure  of  the  dusty-gas  flow  under  supersonic 
conditions . 


Figure  6.  Velocity  and  Density  Distribution  of  a  Dusty-gas 
Shock  Tube 

1.  velocity  distribution 

2.  gas 

3.  density  distribution 

4.  gas 

5.  particle 
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Figure  7.  Variation  of  Pitot  Pressure  with  Time  of  a  Dusty-gas 
Flow 

1.  d  =  10  micron 

Figures  8  and  9  show  the  variation  of  Pitot  pressure  with 
time  at  Mg  =  1.6  and  a  =  0.4  and  1.0.  We  can  see  an  interesting  / 54 
phenomenon,  which  is  that  there  are  many  small  fluctuations 
superpositioned  onto  the  unsteady  part  of  the  Pitot  pressure.  As 
the  cause  was  investigated,  it  was  discovered  that  when  the 
equilibrium  Mach  number  in  the  equilibrium  region  Me<1,  this 
phenomenon  exists.  When  Mg>1,  these  small  fluctuations 
disappear. 


Figures  10  and  11  represent  the  variation  of  the  steady 
value  of  Pitot  pressure  with  <*  at  M  =  1.6  and  3.6,  respectively. 

b 

One  can  see  that  (1)  the  analytical  solution  (equation  (22))  of 
the  effective  flow  in  the  equilibrium  region  is  very  close  to  the 
constant  Pitot  pressure  in  non-equilibrium.  However,  there  is  a 
region  in  Figure  11  where  the  analytical  equation  (equation  (21)) 
is  very  different  from  the  numerical  value.  This  region 
corresponds  to  <*<0.6,  i.e.,  Me<1.  (2)  The  steady  value  of  Pitot 

pressure  increases  with  increasing  <*. 


Figure  10.  Variation  of  Pitot  Pressure  with  <*,  M  =1.6 

s 

1.  numerical  solution 

2.  analytical  solution 
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Figure  11.  Variation  of  Pitot  Pressure  witha,  Mg=3.6 

1 .  numerical  solution 

2.  analytical  solution 

V.  Conclusions 

Through  the  above  discussion,  we  can  reach  the  following 
conclusions.  (1)  Details  of  the  effect  of  the  presence  of  the 
particles  in  a  dusty-gas  shock  tube  on  the  flow  can  be  calculated 
by  using  the  modified  FLIC  method.  This  indicates  that  the 
modified  method  is  reliable.  (2)  The  pressure  of  numerous 
particles  in  the  flow  increases  the  Pitot  pressure.  The  Pitot 
pressure  when  Mg=1.6  and  <*=2  is  1.8  times  that  of  ®=0.  When  M 
=3.6,  it  increases  by  3.1  times.  From  this  we  can  see  that  the 
effect  of  particles  on  Pitot  pressure  also  increases  with 
increasing  Mach  number  Mg .  (3)  When  the  incoming  flow  is  steady, 

Pitot  pressure  reaches  a  steady  value  after  a  relaxation  period. 
When  Mg>1,  the  steady  Pitot  pressure  value  can  be  calculated 
approximately  by  using  the  analytical  solution  (equation  (22))  of 


the  effective  flow  in  equilibrium. 
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The  trial  calculation  is  used  to  explore  some  of  the  features  of  FLIC. 
The  original  approach  to  the  mass  flow  calculation  of  FLIC  is  modified 
through  the  calculation.  The  modified  FLIC  method  is  used  to  compute 
the  dusty-gas  shock  tube  flow  and  its  pitot  pressure.  The  results  obtained 
here  show  clealy  both  the  decay  of  a  frozen  shock  wave  and  a  contact 
surface  and  other  details  of  the  flow.  The  results  also  show  the  variation 
ot  the  pitot  pressure  with  loading  ratio  and  time.  The  stationary  part  of 
the  pitot  pressure  is  in  good  agreement  with  the  analytical  value  of  the 
effective  flow  in  the  equilibrium  region  of  a  dusty-gas  shock  tube  flow. 


An  Experimental  Investigation  of  Flap  Turbulent  Heat  Transfer 
and  Pressure  Characteristics  in  Hypersonic  Flow 
Gao  Ruifeng  (China  Aerodynamic  Research  and  Development  Center) 

Abstract 

This  paper  introduces  the  experimental  results  of  the  heat 

transfer  and  pressure  characteristics  of  a  flap  installed  on  a 

blunt  cone  model  obtained  in  a  shock  tunnel.  Effects  of  flap 

deflection  angle,  attack  angle,  Mach  number  and  unit  Reynolds 

number  on  heat  flow  and  pressure  characteristics  are  discussed. 

The  results  indicate  that  flap  deflection  angle  and  Mach  number 

are  determining  factors  affecting  the  flap  heat  transfer, 

pressure  characteristics  and  control  effectiveness.  The  paper 

also  gives  the  correlation  between  the  peak  heat  flow  and  heat 

pressure,  q  /q  =  (p  / p  )^'^,  as  well  as  an  empirical  formula 

nicix  c  m&x  c 

to  calculate  peak  heating. 

I.  Introduction 

The  three-dimensional  flow  and  its  separation 
characteristics  near  one  flap  on  a  conical  body  are  interesting 
problems  for  us.  When  the  flap  deflects,  it  usually  results  in 
the  interference  between  shock  wave  and  the  boundary  layer  which 
frequently  leads  to  the  formation  of  a  separation  flow. 
Interference  separation  can  change  the  apparent  pressure 
distribution  and  heat  distribution  on  the  aircraft  as  well  as  on 
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the  flap.  Consequently,  it  may  affect  the  stability  and  control 
effectiveness  of  the  aircraft,  as  well  as  may  cause  serious  local 
heating  of  the  surface.  Therefore,  it  is  necessary  to 
investigate  this  type  of  complicated  attached  and  separated 
flows . 

The  re-entry  model  of  a  hypersonic  vehicle  is  shown  in 
photograph  1.  Experimental  studies  on  similar  models  were 

rii 

reported  in  references  [1]  and  [2].  Berman  et  al  used  sharp 
and  little  blunt  cone  models  to  measure  the  heat  transfer, 
pressure  and  drag  increase  on  large  deflection  angles  (45°,  60° 
and  75°)  interference  sweeps  at  M  arl  1  and  Re  =2x10^  (based  on  the 
length  of-  the  model).  Kim  and  Parkinson  used  Calspan's  shock 
tunnel  to  measure  the  heat  transfer,  pressure  and  frictional  drag 
of  the  flap  behind  a  slender  cone  at  M  =  7 . 8''-1 3 . 3  and  Re  =3.3x10~ 
70x10^/foot.  However,  the  flap  deflection  angles  were  small,  10° 
-30°.  The  purpose  of  this  experimental  study  is  to  adopt  this 
model  (see  Photograph  1)  to  increase  the  flap  deflection  angle  to 
measure  the  heat  flow  and  pressure  on  the  flap  and  the  cone 
surface  in  order  to  determine  the  effect  of  factors  such  as  flap 
deflection  angle,  Mach  number,  Reynolds  number  and  attack  angle. 


II.  Experimental  Conditions 
This  experiment  was  conducted  in  a  shock  tunnel.  The 
experimental  gas  is  nitrogen.  The  Mach  numbers  of  the  flow  are 
7.15  and  10.97. 

This  paper  was  received  on  May  16,  1984 
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The  corresponding  Reynolds  numbers  are  4.34x10  m  and  4.14x10  / 5 7 

m  ,  respectively .Experimental  observation  of  a  great  deal  of 
boundary  layer  status  data  indicates  that  the  flap  is  totally  in 
turbulence . 

The  half  conical  angle  of  the  model  is  10°.  The  bluntness 
ratio  is  0.061.  The  win9  deflection  angle  6f  is  22^-51°.  The 

cC 

angle  of  attacks  s  -6°~10°.  The  angle  between  the  flat  surface 
behind  the  conical  model  and  the  axis  of  cone  is  5.5°.  The  two 
flap  chord  lengths  are  80  and  36mm,  respectively.  The  spans  are 
94mm  and  37.5mm.  The  heat  flow  rate  and  pressure  were  measured 
with  a  thin  film  platinum  thermister  and  a  piezoelectric  pressure 
transducer,  respectively. 


III.  Results  and  Discussion 

All  experimental  results  are  expressed  in  terms  of  the  heat 
flow  rate  ratio  (q/q  )  and  pressure  ratio  (p/p  )  where  the 

t  L 

reference  qc  and  pc  are  the  surface  heat  flow  rate  and  pressure 
in  front  of  the  shaved  flat  surface  on  the  rear  of  the  cone  where 
there  is  no  interference. 

The  flow  rate  near  the  flap  of  the  aircraft  is  three- 
dimensional.  When  the  flap  is  deflected  and  the  interference 
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increases,  the  strong  pressure  gradient  causes  the  separation  of 
the  attached  boundary  layer,  the  shock  wave  and  the  interference 
between  boundary  layer  and  the  shock  wave.  Therefore,  the  flow 
field  is  even  more  complicated  (see  photographs  2,  3  and  4).  The 
pressure  distribution  varies  violently  on  and  near  the  flap 
surface  and  the  heating  rate  increases  significantly.  Many 
characteristics  of  this  spatial  flow  are  related  to  various 
influencing  factors  (such  as  flap  deflection  angle,  span,  chord 
length,  flap  profile,  tip  profile,  crevice,  Mach  number,  attack 
angle  and  Reynolds  number).  Effects  of  various  factors  on  the 
thermal  environment  are  described  and  discussed  in  the  following. 
Heat  Transfer  and  Pressure  Characteristics  on  Flap  Surface.  At  M 

00 

=7.15  and  10.97  and  Reao=4.24x10“^m-^  ,  the  flap  heat  transfer  rate 
and  pressure  distribution  at  various  deflection  angles  are  shown 
in  Figure  1.  Apparently,  the  heat  transfer  and  pressure  exhibit 
both  a  boundary  layer  and  a  separation  flow  pattern,  in  agreement 
with  the  general  that  the  flow  is  attached  or  separated  when  the 
interference  is  weak  or  strong.  As  compared  to  the  results 
obtained  with  a  flap  on  a  plate,  the  peak  positions  of  heat 
transfer  and  pressure  are  lagging  behind  for  the  flap  on  a  cone. 
They  vary  between  507*  chord  length  to  807.  chord  length.  When  the 
deflection  angle  is  around  36°,  the  heat  flow  distribution  in  the 
mid  chord  is  slightly  distorted.  This  may  be  due  to  the 
interference  of  the  expansion  wave  at  the  shaved  flat  surface. 

The  following  correlation  formula  is  satisfied  by  the  peak  heat 
flow  and  peak  pressure  in  reference  to  the  rear  cone  values: 
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where  n  =  0.7.  The  correlation  power  index  found  by  Berman  et  al 
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is  n  =  0.8  0.86,  and  n  =  0.8  by  Kim  .  The  difference  may  be 
due  to  the  difference  in  the  local  geometrical  profile  of  the 
model.  Therefore,  we  believe  that  the  correlation  index  is 
determined  by  the  local  flow  field  characteristics. 
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Figure  1.  Heat  Flow  Distribution  and  Pressure  Distribution  at 
the  Center  of  the  Flap  Wing 

1.  Re®  =  4.34  x  lO^rn”! 

2.  Re®  =  4.14  x  lum” 

3.  heat  flow  distribution 

4.  Re®  =  4.34  x  104m 

5.  Re®  =  4.14  x  10'm 

6.  pressure  distribution 

From  an  engineering  angle,  on  the  basis  of  the  analysis 
given  in  reference  [2]  and  by  using  the  experimental  results  in 
this  work,  the  following  empirical  formula  to  calculate  the  peak 
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heat  transfer  is  obtained 


L=s  =  9.136-18.39  <3,- (23.195  -  56.08  d,> 
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where  6^  is  in  radians  and  Re^  is  the  unit  Reynolds  number  (m~‘).  / 58 

The  applicable  range  is:  28°<6£^45°,  6c=8 ,5°^1 1  °  ;  Mob=7^1  1 ;  and 
small  bluntness.  The  error  in  most  cases  is  less  than  +10X.  The 
formula  shows  that  the  peak  heat  flow  is  proportional  to  the  flap 
deflection  angle  and  Mach  number,  and  inversely  proportional  to 
Re 

<3D 

Characteristics  of  Heat  Flow  in  Separation  Region.  The 
deflection  of  the  flap  causes  the  thermal  environment  near  the 
flap  to  vary  significantly.  Observation  of  the  shadow  (see 
photographs  2,  3  and  4)  and  measurement  of  the  heat  flow  indicate 
that  the  separation  shock  wave  shifts  forward  with  increasing 
deflection  angle.  The  region  influenced  by  separation  is 
expanding  (see  Figure  2).  v 


,  33*48'  39*18'  50*3077  . 

i,t  •  *  •  ■  •  *  f!  p  _S_ 
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Figure  2.  Heat  Flow  Dis tiioution  in  Separation  Region 

(Mm  *7.15  .  Rrm  =  4.34x  10’  /|| 
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(photographs) 

Effect  of  Attack  Angle:  Under  a  specific  separation 

condition  (6^=39°16,)>  the  attack  angle  of  the  model  was  changed 

to  -4°,0°,  4  and  10°  to  measure  the  heat  flow  and  pressure  on  the 

flap.  The  results  are  shown  in  Figure  3.  It  was  found  that  the 
within  a  small  angle  of  attack  range 

heat  flow  varies  slightly^  <-4°<a  <4")t  Hence,  the  heat  transfer 
ratio  at  zero  attack  angle  can  be  roughly  used  to  represent  the 
heat  transfer  ratio  in  the  range  of  small  attack  angles. 

However,  when  the  attack  angle  is  over  10°,  the  heat  flow  on  the 
flap  is  decreased  significantly. 
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Figure  3.  Effect  of  Attack  Angle  on  Flap  Heat  Flow  Distribution 

Under  the  condition  of  attached  flow  (6£=22.5°),  the  heat 

transfer  ratio  and  pressure  ratio  on  the  flap  decrease  as  the 

attack  angle  decreases  from  0°  to  -4°  and  -6°,  i.e.,  the  flap  is 

the  leeward  surface.  The  flap  heat  and  pressure  increase  must 

be  slower  than  the  corresponding  values  on  the  rear  cone  surface. 

Effect  of  Reynolds  Number.  In  this  work,  the  Mach  number 

was  kept  unchanged,  i.e.,  M^IO.84.  When  the  unit  Reynolds  / 59 

7-1  7-1 

number  Rea#=2. 82x10  m  and  4.14x10  m-  ,  the  measured  flap  heat 
transfer  ratio  and  pressure  ratio  basically  do  not  vary. 

Therefore,  we  believe  that  the  dependence  of  characteristics 
controlling  the  flap  heat  transfer  ratio  and  pressure  ratio  on 
unit  Reynolds  number  is  weak  in  the  experimental  conditions 
studied  (see  Figure  4). 

Effect  of  Mach  Number.  When  the  flap  deflection  angle  is 

7  - 1 

36.5°  and  the  Reynolds  number  is  approximately  4.24x10  m”  ,  at 
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10.97  and  7.15,  experimental  results  show  that  the  heat 
transfer  ratio  and  pressure  ratio  on  the  flap  are  strongly 
dependent  on  the  Mach  number.  With  increasing  Mach  number,  the 
heat  transfer  and  pressure  ratio  will  increase  (see  Figure  5). 


Figure  4.  Effect  of  Unit  Reynolds  Number  on  Flap  Heat  Flow  and 
Pressure  7  - 

/"4.14  x  10m-  (hollowlN 

^<■1 0.84, Re  *  1  7  1  / 

V2.82  x  10m-' (solid)  J 
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1.  heat  flow 

2.  pressure 


Figure  5.  Effect  of  Mach  Number  on  Flap  Heat  and  Pressure 
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1 .  heat  flow 

2.  pressure 


IV.  Conclusions 

Under  the  experimental  flow  condition,  we  reach  the 
following  conclusions  in  studying  the  turbulent  thermal 
environment  on  a  flap  in  a  blunt  cone  model: 

1.  The  heat  flow  and  pressure  are  distributed  in  a  boundary 


layer  or  separation  flow  pattern  on  the  flap.  The  peak  heat 
transfer  and  peak  pressure  (in  reference  to  the  value  on  the  rear 
body)  satisfy  the  following  correlation: 


2.  The  variation  of  flap  re-attached  flow  with  attack  angle 
is  similar  to  that  of  rear  cone  heat  flow  with  attack  angle  (in 
the  small  attack  angle  region). 

3.  The  re-attached  heat  flow  and  pressure  characteristics 
of  the  flap  are  strongly  dependent  on  the  Mach  number.  Their 
dependence  on  the  unit  Reynolds  number  is  weak. 
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AN  EXPERIMENTAL  INVESTIGATION  OF  FLAP  TURBULENT 
HEAT  TRANSFER  AND  PRESSURE  CHARACTERISTIC'S 
IN  HYPERSONIC  FLOW 

Gao  Ruifeag 

(China  Aerodynamic  Research  and  Deveiorment  Center) 

Abstract 

This  paper  presents  experimental  results  of  flap  heat  transfer  and 
pressure  characteristics  on  a  blunt  cone  in  shock  tunnel.  Effects  of 
flap  deflection  angle,  angle  of  attack,  Mach  number  and  unit  Reynolds 
number  are  discussed. 

Results  have  shown  that  flap  deflection  angle  and  Mach  number  are 
decisive  factors,  which  considerably  affect  the  flap  heating,  pressure 
characteristics  and  control  effectiveness.  This  paper  gives  a  correlation 

of  peak  heating  and  peak  'pressure,  j  <  also  gi'«s  an  empi- 

rical  folmula  for  estimating  peak  heating. 
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LDA  Measurements  for  Leading  Edge  Vortex  Core  Velocity  of  a 

Strake-wing 

Lu  Zhiyong  (Beijing  Institute  of  Aeronautics  and  Astronautics) 
and  Cheng  Yuangzhong  (Institute  No.  304  of  Aviation  Ministry) 

Abstract 

The  axial  velocity  distribution  of  a  vortex  core  on  a 
strake-wing  was  measured  by  a  two-dimensional  laser  Doppler 
anemometer  (LDA)  in  a  water  funnel.  Furthermore,  the  spatial 
position  of  the  strake-wing  vortex  and  the  variation  of  the 
vortex  breakdown  point  with  the  angle  of  attack  are  given. 

It  is  demonstrated  in  this  paper  that  the  axial  velocity 
distribution  at  various  attach  angles  can  be  expressed  by  one 
curve  by  choosing  an  appropriate  non-dimensional  parameter.  But, 
this  is  not  true  at  downstream  of  the  vortex  breakdown  point 
because  of  the  violent  change  in  velocity. 

Introduction 

Over  a  long  period  of  time,  the  aerodynamic  characteristics 

and  the  flow  phenomena  of  slender  delta  wings  and  strake-wings 

-ri  o  31 

have  been  widely  studied  ’  It  was  found  that  the 

aerodynamic  characteristics  of  these  wings  at  large  angles  of 
attack  are  determined  by  the  leading  edge  vortex  (see  Figure  1). 
As  we  all  know  very  well  that  the  leading  edge  vortex  originates 
from  the  separation  of  the  shear  layer  at  the  leading  edge  of  the 


wing.  Most  of  the  vortex  moment  is  concentrated  in  the  vortex 
core.  The  flow  field  of  a  vortex  core  has  been  studied 
extensively.  However,  most  of  the  studies  limited  to  the 
measurement  of  heat  and  direction^’^  .  Different  from  those 
methods  mentioned  before,  LDA  shows  its  unique  advantage  that  it 
will  not  interfere  with  the  flow  field^^.  In  addition,  its 
measuring  signal  is  transmitted  at  the  speed  of  light.  In 
conjunction  with  a  processing  system  including  a  computer,  it  is 
capable  to  perform  real  time  measurement.  This  is  very  useful  to 
the  study  of  a  vortex  core  flow  field  and  the  measurement  of 
turbulent  parameters  in  the  vortex  core.  The  accurate 
measurement  of  the  flow  field  can  further  deepen  the  study  on  the 
leading  edge  vortex,  which  will  facilitate  the  creation  of  more 
appropriate  mathematical  models. 


Figure  1.  Leading  Edge  Vortex  of  a  Delta  Wing 

Experimental  Apparatus  and  Model 

Water  Tunnel 

The  measurement  of  the  flow  field  of  the  vortex  core  of  a 
strake-wing  was  performed  in  the  water  tunnel  at  Beijing 
Institute  of  Aeronautics  and  Astronautics.  The  experimental 
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section  of  the  water  tunnel  is  0.4x0. 4m  in  square  cross-section. 
It  is  6m  long.  The  Reynolds  number  based  on  the  average  chord 
length  of  the  model  as  the  reference  is  Re=12000. 

This  paper  was  received  on  June  10,  1984 
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Hydrogen  Bubble  Generator 

The  principle  is  shown  in  Figure  2.  When  the  aluminum  model 
is  placed  between  the  platinum  wire  and  the  carbon  rod,  the  model 
itself  becomes  an  electrode.  Hydrogen  bubbles  are  generated  at 
its  leading  edge.  When  the  model  is  at  a  certain  angle  of 
attack,  the  hydrogen  bubbles  display  the  core  of  the  leading  edge 
vortex.  When  we  perform  the  LDA  measurement,  hydrogen  bubbles 
are  used  as  the  tracing  particles.  The  size  of  the  hydrogen 
bubble  is  usually  below  4  microns. 


Figure  2.  Hydrogen  Bubble  Generator 


1 .  input  end 

2.  transformer 

3.  relay 

4.  to  metal  wire 

5.  polarity  switch 


The  DISA  two-dimensional  Laser  Doppler  Anemometer  (LDA) 

The  LDA  system  is  shown  in  Figure  3.  In  the  experiment,  the 
optical  system  of  the  LDA  is  shown  in  Figure  4. 
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Figure  3.  Block  Diagram  of  the  LDA  System 

1 .  laser 

2.  beam  modulator  and  converter 

3.  light  detector 

4.  signal  processing  display 

5.  receiving  PDP11/03 

6.  terminal 
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Figure  4.  Optical  Arrangement  of  the  LDA  System 


2.  blue 

3.  green 

4.  mirror  M 

5.  polychromatic  light 

6.  blue  and  green  light 

7.  optical  element 

8.  laser 

9.  cover  plate  P 

10.  receiver 

Model 

The  profile  and  size  of  the  strake-wing  is  shown  in  Figure 
5.  The  sweepback  angle  of  the  leading  edge  of  the  strake  is  X=76 
.  The  sweepback  angle  of  the  primary  wing  is  X=30°. 


Figure  5.  Profile  of  the  Strake-wing 

Results  and  Discussions 

1.  Position  of  Vortex  Core  and  Breakdown  of  Vortex  Core 
Figure  6  shows  the  positions  of  the  vortex  core  at  various 
angles  of  attack,  which  were  obtained  through  the  use  of  hydrogen 


bubbles  and  laser  measurements.  Figure  7  shows  a  comparison  of 
the  breakdown  point  of  the  vortex  core  for  a  strake-wing  to  that 
of  a  delta  wing.  In  the  current  experiment,  the  breakdown  point 
has  already  reached  the  trailing  edge  when  the  attack  angle  a=18° 
This  is  slightly  lower  than  the  result  obtained  from  wind  tunnel 
measurement.  It  should  be  pointed  out  that  the  position  of  the 
breakdown  point  does  not  remain  unchanged  after  it  is  moved  onto 
the  primary  wing.  On  the  contrary,  the  breakdown  points  of  the 
two  primary  vortices  are  located  on  either  side.  Each  of  them  is 
mo\ing  back  and  forth  around  a  specific  average  position.  The 
phase  of  this  movement  differs  by  180°. 

'  2.  Axial  Velocity  Distribution  of  the  Vortex 
According  to  reference  [8],  the  space  occupied  by  the 
leading  edge  vortex  can  be  divided  into  three  regions,  as  shown 
in  Figure  8.  The  parameters  D  and  \/ vt  are  the  major 
characteristics  of  these  three  regions.  D  represents  the 
distance  between  the  spiral  vortex  sheets ,  v' vt  Is  the  viscous 
diffusion  distance  and  v  is  the  dynamic  viscosity  coefficient,  t 
is  the  characteristic  time  which  can  be  expressed  as  x/U^,  where 
x  is  the  coordinate  in  the  direction  of  the  primary  flow  and 
is  the  incoming  flow  velocity.  The  viscous  region  is  D<< VrvF. 

In  this  region,  the  vortex  core  (also  called  the  subcore)  is  also 
spinning  as  a  solid.  There  is  only  axial  velocity  and  no 
rotational  velocity  on  the  axis  of  the  vortex. 


Figure  6.  Position  of  Vortex  for  a  Strake-wing 


1.  Strake-wing  Tunnel  Experiment 

2.  75°  Sweepback  Delta  Wing  Water  Tunnel  Experiment 

3.  Hydrogen  Bubble  in  Water  Tunnel,  Strake-wing 

Figure  9  shows  the  axial  velocity  distribution  of  the  vortex 
core  at  4  different  angles  of  attack.  When  a=12°,  based. on 
hydrogen  bubble  observation,  the  breakdown  point  is  far  away  from 
the  trailing  edge.  The  axial  velocity  begins  to  increase  from 


the  apex.  This  is  because  that  vortex  moment  is  continuously 
added  to  the  vortex  core  along  the  leading  edge  from  the  top. 
Their  trajectory  can  be  considered  as  a  series  of  spiral  lines. 

An  axial  velocity  component  pointed  downstream  is  created  with 
respect  to  the  vortex  center.  Near  the  trailing  edge,  because  of 
the  inverse  pressure  gradient,  the  axial  velocity  of  the  vortex 
core  drops.  When  a=18°,  the  breakdown  point  of  the  vortex  core 
has  already  reached  the  trailing  edge.  In  this  case,  the  axial 
velocity  of  the  vortex  center  drops  abruptly  behind  the  breakdown 
point.  The  deceleration  gradient  is  very  large.  It  is  not  the 
same  as  the  deceleration  gradient  near  the  trailing  edge  at  a=12° 
.  At  a=25°,  the  breakdown  point  has  already  moved  forward  to  the 
wing  surface,  approximately  at  817.  of  the  root  chord.  At  a=35°, 
the  vortex  breakdown  point  has  moved  to  52%  of  the  root  chord. 
Because  the  vortex  intensity  increases  and  the  induction  effect 
of  the  vortex  intensifies,  the  negative  pressure  on  the  wing 
(relative  to  the  static  pressure  of  the  incoming  flow)  increases, 
and  the  maximum  axial  velocity  also  increases.  The  axial 
velocity  can  be  2.4  times  the  incoming  flow  velocity.  It  should 
be  pointed  out  that  this  value  is  lower  than  that  for  a 
corresponding  delta  wing.  This  is  because  the  primary  wing  of 
the  strake-wing  cannot  further  reinforce  the  vortex. 


Figure  8.  Structure  or  the  Vortex 


Figure  9.  Axial  Velocity  Distribution  Along  the  Center  of  the 
Vortex 

3.  Similarity  Curve  of  Axial  Velocity  Distribution  at  the  /64 
Vortex  Center 

Upstream  from  the  vortex  breakdown  point,  the  curves  of 
axial  velocity  vs.  attack  angle  a  have  some  similarities.  At  a 


fixed  angle  of  attack,  the  maximum  axial  velocity  Uc  is  used  as 
the  reference  velocity  to  render  velocity  non-dimensional.  The 
distance  between  the  maximum  velocity  point  on  the  vortex  axis 
and  the  apex  of  the  wing,  Sm>  is  used  as  the  family  of  curves 
shown  in  Figure  9  can  be  expressed  as  the  curve  in  Figure  10. 
However,  it  is  not  possible  to  express  it  downstream  from  the 
vortex  breakdown  point  by  one  curve.  Whether  this  property  has  a 
more  generalized  significance  remains  to  be  further  verified. 


4.  Velocity  Variation  at  the  Vortex  Breakdown  Point 
The  breakdown  phenomenon  itself  is  also  on  unsteady  process. 
As  described  before,  the  position  of  the  breakdown  point  varies 
with  time.  The  velocity  near  the  breakdown  point  also  varies 
randomly.  Table  1  shows  the  mean  velocity  distribution  in  a  very 
short  period  of  time  at  the  breakdown  point.  It  also  shows  the 
corresponding  mean  square  root  velocity  values.  Therefore,  we 
should  measure  dynamic  parameters  in  order  to  understand  the 
breakdown  phenomenon  well  and  to  know  the  dynamic  characteristics 


of  the  flow  field  after  the  breakdown. 


For  delta  and  strake-wings  with  a  sharp  leading  edge,  the 
major  parameters  affecting  the  breakdown  characteristics  are  the 
angle  of  attack  a  and  the  sweepback  angle  X.  The  effect  of 

[Q] 

Reynolds  number  is  secondary  .  The  direct  cause  for  vortex 
breakdown  is  the  effect  of  viscosity  and  inverse  gradient  in  the 
vortex  core.  This  is  a  commonly  accepted  viewpoint. 


Table  1  Velocity  and  Its  Mean  Square  Root 
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1.  sampling  time  sequence* 

2.  average  velocity  within  0.1  sec,'UA/U 

3.  root  mean  square  within  0.1  sec,  Uv/uX 

4.  *  The  time  in  each  time  sequence  number  is 
0.1  sec.  There  are  100  samplings 


Conclusions 

1.  With  increasing  attack  angle,  the  angle  between  position 
of  the  vortex  axis  and  the  wing  surface  increases  and  the  angle 
between  the  plane  of  symmetry  of  the  wing  and  the  vortex  axis 
decreases.  This  is  an  important  factor  for  the  presence  of  the 
non-linear  (upward)  torque  on  a  single  wing. 

2.  The  axial  velocity  on  any  point  on  the  vortex  axis 
begins  from  the  apex  and  reaches  its  maximum  at  a  point  in  front 


of  the  trailing  edge.  Then,  it  begins  to  decrease.  At  the  wing 
surface  where  breakdown  takes  place,  the  decrease  of  velocity  is 
caused  by  an  inverse  pressure  gradient,  leading  to  the  breakdown 
of  the  vortex. 

3.  Appropriate  non-dimensional  parameters  are  used  to 
express  the  velocity  distribution  at  any  point  along  the  vortex 
axis  of  a  strake -wing  and  can  be  expressed  as  one  curve.  Beyond  the  /65 
breakdown  point,  this  general  curve  does  not  exist. 

In  the  research  process,  some  comrades  from  the  506  teaching 
and  research  offices  and  the  304  Institutes  assisted  in  the  work. 
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LDA  MEASUREMENTS  FOR  LEADING  EDGE 
VORTEX  CORE  OF  A  STRAKE-WING 

LQ  Zhiyong 

(Billing  Institute  of  Aeronautics  eng  Attronau tict) 

Cheng  Yuanzbong 
(Institute  No.  304  of  Aviation  Ministry) 

Abstract 

-  X 

The  axial  velocity  distribution  along  leading  edge  vortex  core  of  a 
strake-wiag  were  measured  by  means  of  2-Dimensional  L.  D.  A  and  hy¬ 
drogen  babble  technique  in  the  water  tunnel.  The  trajectories  of  the  vortex 
core  aad  the  vortex  breakdown  point  which  vary  with  angles  of  attack 
are  preseated.  It  has  been  showa  that  axial  velocity  distribution  at  dif¬ 
ferent  angles  of  attack  can  be  fitted  with  a  aondimensional  curve.  But  this 
is  not  Araf  at  dowasteam  of  vortex  breakdown  point  because  of  dramatic 
chaagf  jot.  the  velocity.  . 
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An  Experimental  Research  of  Boundary  Layer  Control  Technique 
for  Low  Speed  Multi-component  Airfoils 
Wang  Jie  (Northwestern  Polytechnical  University) 

I.  Introduction 

In  order  to  improve  the  takeoff  and  landing  characteristics 
of  an  aircraft,  it  is  necessary  to  conduct  low  speed  multi- 
component  wind  tunnel  experiments  in  the  design  and  selection  of 
various  devices  to  enhance  lift. 

The  height  to  width  ratio  of  a  conventional  two-dimensional 
experimental  section  is  relatively  large  and  the  geometric  aspect 
ratio  of  the  model  is  small.  The  boundary  layer  on  the  wall  in 
the  experimental  section  will  affect  the  experimental  result. 

The  narrower  the  model  is,  the  more  serious  the  effect  becomes. 
The  turbulent  boundary  layer  of  the  tunnel  wall  expands  toward 
the  laminar  flow  region  at  the  leading  edge  of  the  model  in  a 
wedge  area,  which  moves  the  turning  point  on  the  middle  cross- 
section  of  the  model  forward.  When  a  pressure  gradient  exists  in 
a  two-dimensional  airfoil  experiment,  the  flow  at  the  wall 
becomes  extremely  complicated.  The  total  intensity  is  equivalent 
to  the  scatter  of  the  vortex  circulation  of  the  airfoil  on  the 
tunnel  wall  in  all  directions.  Induced  velocity  and  induced 
attack  angle  are  created  on  all  cross-sections  in  the  airfoil. 
Consequently,  the  two-dimensional  nature  of  the  flow  is  affected. 

The  most  important  effect  is  the  premature  separation  at  the 
model-tunnel  wall  junction.  Far  before  separation  takes  place  on 
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the  middle  cross-section  of  the  model,  the  boundary  layer  near 
the  model-tunnel  wall  junction  is  separated  under  a  relatively 
small  inverse  preessure  gradient.  Furthermore,  it  propagates 
approximately  along  a  45°  angle  toward  the  middle  of  the  model. 
For  models  with  an  aspect  ratio  less  than  2,  premature  flow 
separation  at  the  tunnel  wall  has  already  affected  the  pressure 

[3]  [2] 

measurement  at  the  middle  cross-section  ’ 

II.  Experimental  Apparatus,  Blowing  System,  Model  and  Experi¬ 
mental  Technique 

The  experiment  was  performed  in  a  direct  flow,  closed,  low 

speed,  two-dimensional  wind  tunnel.  The  cross-section  in  the 

experimental  section  is  rectangular.  It  is  2m  high  and  0.2m 

wide.  The  boundary  layer  thickness  on  the  side  wall  of  the 

r  3 1 

experimental  section  was  measured  to  be  30mm  .  The 
experimental  wind  speed  is  V=45m/sec  and  the  Reynolds  number 
Re=1.2x106. 

The  blowing  system  has  6  blowing  slots  which  are  installed 
symmetrically  on  both  sides  of  the  turning  wheel  (See  Figure  1). 
The  front  slot  is  located  on  the  upper  surface  of  the  leading 
edge  of  the  airfoil.  It  is  150mm  long.  The  middle  slot  is 
located  at  65%  chord  length  from  the  leading  edge.  It  is 
perpendicular  to  the  wing  chord  and  symmetric.  The  total  length 
is  330mm.  The  rear  slot  is  located  at  the  upper  surface  of  the 
slap.  It  is  80mm  long. 

Each  blowing  slot  is  2mm  wide.  Compressed  air  is  blowing 
out  of  a  pressure  stabilizer  along  the  tunnel  wall  surface.  It 


was  measured  that  more  than  907«  of  the  flow  blowing  out  of  the 
slot  along  the  blowing  direction  is  homogeneous.  Furthermore, 
the  flow  blowing  from  the  slots  to  eliminate  the  wall  boundary 


layer  does  not  have  a  significant  effect  on  the  tunnel  flow  field 
itself C3] . 

This  paper  was  received  on  December  27,  1983  and  revised  manu¬ 
script  was  received  on  April  16,  1984. 

The  experimental  models  are  single  and  double  slotted  flap 
airfoils.  There  are  54  pressure  measuring  holes  on  the  model. 

The  chord  length  of  the  basic  airfoil  is  400mm. 

The  blowing  amount  is  controlled  by  a  blowing  coefficient 
which  is  expressed  as: 

7-  =  <?  (1) 

Q=Po-P.  (2) 

7.-C p.,-p.)!qm  (3) 

where  — dynamic  pressure  at  the  outlet  of  the  blowing  slot,  qo- 

-dynamic  pressure  in  the  experimental  section;  and  po — static 
pressure  in  the  experimental  section.  There  is  a  pressure 
measuring  hole  on  the  side  of  the  pressure  stabilizing  box.  The 
pressure  sensed  through  the  hole  is  considered  to  be  the  total 
pressure  of  the  blowing  jet  PQj •  This  pressure  is  connected  to  a 
pressure  gauge.  The  pressure  gauge  reaching  is  used  to 
accurately  control  the  blowing  quantity. 


Figure  1.  Blowing  Slot  Installed  on  the  Upper  Surface  of  the 
Model 

According  to  the  measured  pressure  distribution  on  the 
surface  of  the  airfoil,  we  can  calculate  the  pressure  coefficient 
Cp.  Then,  through  numerical  integration  and  tunnel  wall 
interference  correction,  it  is  possible  to  obtain  the  lift,  drag 
and  torque  coefficient  of  the  whole  airfoil.  The  airfoil  drag  is 
obtained  by  measuring  the  total  pressure  loss  in  the  wake. 

III.  Experimental  Results  and  Discussion 

From  the  part  of  experimental  curve  already  obtained  (See 

Figure  2)  we  can  see  that  the  slope  of  the  lift  curve  and  the 

maximum  lift  coefficient  C  are  relatively  low  when  not 

ymax  3 

blowing.  When  the  wall  boundary  layer  was  blown  away  by  blowing 
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an  amount  of  air  corresponding  to  An  and  Ag  equal  to  10,  Cy  and  Cymax 
are  apparently  larger.  This  blowing  quantity  already  made  the 
flow  separation  appear  in  the  middle  cross-section  of  the 
model,  instead  of  at  the  junction  of  the  mode  and  the  tunnel 
wall.  Blowing  through  the  front  slot  is  more  effective  in 
preventing  separation  of  .the  tunnel  wall  boundary  layer.  When 
the  blowing  coefficient  exceeds  15,  Cy  essentially  remains 
unchanged  as  the  blowing  quantity  increases.  In  other  words, 
once  the  blowing  exceeds  the  required  quantity  to  avoid  the 
premature  separation  at  the  model-tunnel  wall  junction,  the 
effect  of  blowing  is  no  longer  obvious.  Based  on  the 
experimen-tal  results  obtained  with  multi-component  airfoils  with 
three  relative  thicknesses  at  16.57.-187.  in  this  work,  the 
appropriate  blowing  coefficient  should  be  controlled  between  15- 
20.  Blowing  away  the  tunnel  wall  boundary  layer  restores  and 
Cymax*  They  are  increased  by  157.-307.  and  0.3-0. 4,  respectively, 
as  compared  to  those  obtained  without  blowing.  ' 


Figure  2.  Effect  of  Blowing  on  the  C^-a  Curve 

X  -front  slot  blowing  coefficient 
n 

x  -middle  slot  blowing  coefficient 
s 

1.  double-slotted  flap  airfoil,  flap  deflection 
angle  <t>  =  0°,  15°,  20° 

Because  the  width  of  the  experimental  section  is  very 
narrow,  when  the  blowing  coefficient  is  too  large,  an  overblowing 
effect  will  appear  in  the  experiment,  which  will  directly  affect 
the  pressure  measurement  at  the  middle  cross-section  of  the 
model.  This  is  reflected  on  the  curve.  The  slope  of  the  lift 


curve  remains  unchanged.  The  curve  is  shifted  upward  and  the 
zero  lift  attack  angle  increases.  Therefore,  the  control  of  the 
blowing  quantity  in  a  small  wind  tunnel  will  affect  the 
reliability  of  the  experimental  result. 

When  comparing  the  pressure  distribution  with  that  without  /68 
blowing, its  slope  remains  unchanged.  The  surface  area  beneath 
the  airfoil  is  essentially  unchanged,  while  the  upper  surface 
area  is  increased  (premature  separation  at  the  tunnel  wall 
affecting  a  larger  region). 

In  an  actual  two-dimensional  flow,  the  pressure  differential 
drag  obtained  from  the  pressure  distribution  should  be  less  than 
the -airfoil  drag  measured  in  the  wake.  The  latter  is  the  sum  of 
pressure  differential  drag  and  friction  drag.  However,  the 
airfoil  drag  when  there  is  no  blowing  is  smaller  than  the 
pressure  differential  drag  (See  Figure  3).  As  described  before, 
the  vortex  from  the  tunnel  wall  changes  the  attach  angle  at  every 
cross-section.  The  induced  drag  is  already  included  in  the 
pressure  measured  at  the  cross-section  in  the  middle  (See 
Figure  3). 


<•>  %n»n6*o'  <*»  = 


Figure  3. 


Effect  of  Blowing  on  Pressure  Differential  Drag  and 
Airfoil  Drag 


(a)  double-slotted  flap  6=0° 

(b)  double-slotted  flap, 6=38' 

1.  differential  drag 

2.  cd  airfoil  drag 


IV.  Conclusion  and  Prospect 


1.  Using  a  tunnel  wall  boundary  layer  control  technique,  i 
is  possible  to  establish  an  approximate  two-dimensional  flow 


state  until  stalling.  Reliable  value  can  be  obtained. 


2.  The  location  of  the  blowing  slot  is  not  the  key  to 
preventing  the  large  area  flow  separation  beyond  the  slot.  (In 
this  experiment,  the  position  of  the  front  slot  was  changed, 
however,  we  could  not  see  any  obvious  effect  on  the  experimental 
result).  Nevertheless,  blowing  slots  perpendicular  to  the  local 


flow  are  optimal. 

3.  Blowing  from  the  slot  is  most  effective  on  suppressing 

the  effect  of  the  tunnel  wall  boundary  layer.  Only  the  front  and 

middle  slots  are  needed  for  satisfactory  experimental  results. 

In  the  double-slotted  flap,  the  blowing  quantity  from  the  front 

slot  should  be  higher  than  that  from  the  middle  slot. 

The  experimental  research  of  using  a  blowing  control 

technique  for  the  boundary  layer  is  still  preliminary.  These 

preliminary  conclusions  are  obtained  based  on  the  experimental 

study  of  airfoils  with  large  relative  thickness.  For  airfoils 

with  a  smaller  relative  thickness  such  as  modern  airfoils 

including-  the  supercritical  airfoils,  what  combination  of  blowing 

quantities  can  be  used  to  obtain  reliable  two-dimensional 

experimental  data  and  whether  a  reliable  C  value  can  be 

ymax 

obtained  by  blowing  and  the  appropriate  correction  of  the  Re 
number  still  remain  to  be  studied  and  explored. 

This  work  received  the  guidance  and  assistance  from  comrades 
Xia  Yushun,  Xi  Zhongxiang  and  Bao  Guhua.  Comrades  in  the 
laboratory  also  actively  coordinated  the  experimental  work.  The 
author  wishes  to  express  his  gratitude  to  them. 
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'AN  EXPERIMENT  RESEARCH  OF  BOUNDARY 
LAYER  CONTROL  TECHNIQUE  FOR 
MULTI-COMPONENT  AIRFOILS 


Wang  Jie 

•V-  •-  ,  iNorthwttUrm  Poly  technical  University) 

*  Abstract 

jr  Jhe  present  paper  describes  some  research  results  of  boundary  layer 
control  on  the  wall  of  a  two  dimensional  low  speed  wind  tunnel  with  blow¬ 
ing  system.  In  the  experimental  research  the  airfoil  models  with  single 


and  double  slotted  trailing  edge  flap  were  adopted.  The  boundary  layer 
control  is  a  very  effective  method  to  prevent  the  premature  flow  sepa¬ 
ration  in  the  junction  corner  between  the  tunnel  wall  and  the  airfoil  model 
and  to  get  rid  of  three  dimensional  effect  for  single  or  multi-component 
airfoil  research. 


In  the  paper,  some  of  contrast  research  results  are  given  by  figures 
and  the  selection  of  blowing  quantity  and  the  disposition  of  the  blowing 
slots  are  described. 
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A  Simpler  Implicit  Method  to  Solve  N-S  Equation 
Zhang  Hanxin,  Yu  Zechu  and  Lu  Linsheng 
(All  of  China  Aerodynamic  Research  and  Development  Center) 


/  70 


1.  Introduction 

In  the  past  decade,  a  great  deal  of  progress  was  made  in 

solving  the  N-S  equation  of  a  separation  flow  by  numerical 

methods.  Explicit  finite  difference  methods  have  been  developed 

to  a  mature  stage.  However,  because  of  the  limitation  of 

stability  conditions,  it  takes  longer  to  calculate  a  stable 

solution.  In  order  to  overcome  this  difficulty,  implicit 

methods,  such  as  the  implicit  factoring  method  and  implicit 

alternating  direction  method,  were  developed.  However,  these 

methods  require  solving  three  diagonal  matrices.  As  compared  to 

explicit  methods,  the  program  for  an  implicit  method  is  more 

complex  and  the  computation  time  needed  for  every  step  is  longer. 

In  order  to  avoid  the  difficulty  in  the  iteration  of  the  three 

r  11 

diagonal  matrices,  MacCormack  introduced  a  new  implicit  method 
and  its  iteration  process  involves  the  chasing  of  two  diagonal 
matrices.  On  the  basis  of  reference  [1],  an  even  simpler 
implicit  method  is  introduced  in  this  paper.  Its  iteration 
process  only  involves  scalar  chasing.  Examples  given  in  this 
work  show  that  this  method  is  capable  of  giving  results  in 
agreement  with  those  obtained  by  explicit  methods.  Its 
efficiency  is  comparable  to  that  of  reference  [1],  but  it  is 
simple  and  easy. 
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2.  Difference  Scheme 
1.  Model  Equation 

Let  us  study  the  following  equation: 

du  ,  dF,  dF,  _Q 

d/  a*  d*  (2.1) 

where 

F,  **•«  (o*0)  (2.2) 

(2.3) 

a  and  v  are  constants.  Let  us  use  the  following  second  order 
accuracy  implicit  difference  scheme: 


-.p  =  uf*T  - C  (1 -r  -  Fif)  -  ^  (F!^3 ,  -  F -f)  ] - 

At 


(2.4) 


Li*  /]?»■»  i  ri-r  u 


-±-  (i»*+«r‘) 

Z 


This  paper  was  received  on  May  14,  1984.  The  content  of  this 
paper  was  introduced  in  the  Fluid  Mechanics  Computation  Class  at 
Qinghua  University  in  June  1983. 

where  3  is  constant  determined  by  o  and  v.  If  the  first  step  is  / 7 1 
to  pick  a  backward  difference,  then  choose  3  so  that  B(dF/9u)  =3a 
>0.  If  the  first  step  is  to  pick  a  forward  difference,  then 
choose  3  so  that  S( 9F/ 9u)=3a<0 .  Let 

Suy7*  ur~‘ 

d«frr=urT 

.F-F.  +  F, 


-V, 


-u, 


(2.5) 


Using  Taylor  expansion,  equation  (2.4)  can  be  written  as: 


VM1 


m 


J* •». 


( ,  ^  ^ )  a.  :v'  -  -  f  ; f  :TI  >  +  <<  •  fi  ^ 


(2.6) 


(duy7+duy‘) 


Here,  d=Ba>0.  If  a  stability  analysis  is  performed  on  equation 
(2.6),  it  is  found  that  in  order  to  satisfy  the  stability 


condition  d  should  be: 


{[-r  ("£)-&]•  •} 


(2.7) 


where  k  is  an  insurance  index  which  is  larger  than  1 .  One  can 
see  that  the  difference  formula  (2.6)  thus  derived  is 
MacCormack's  implicit  formula.  When  d=0 ,  it  is  MacCormack's 
explicit  formula. 

2.  One-dimensional  Unsteady  N-S  Equation 


This  equation  is 


+  (2.8) 

dt  dx  ox 


Here,  U  is  an  n  dimensional  vector,  F<j=F^(u)  is*  a  function  of  a 
n-dimensional  vector  which  is  not  dependent  on  viscosity,  and  F 
is  a  function  of  an  n-dimensional  vector  which  is  dependent  on 
viscosity.  Obviously,  equation  (2.8)  can  be  written  as: 


dU  a  dU  +  dF±  =  o 
dt  dx  dx 


(2.9) 


where  A  =  (aF^/aU).  Let  us  assume  that  it  is  possible  to  express 

- 1  11  - 1 

A  as  A=SaS  .  Here,  A=diag(x,)  and  is  the  eigenvalue  of  A.  S 

A  A 

is  the  lift  eigen  vector  matrix  of  A  and  S  is  its  inverse  matrix. 
In  order  to  solve  equation  (2.9),  in  analogy  to  the  model 


**<•*»  *  *  •  ■  a  '  m  *  ■  f  »  ™  i»  *  **»*►*-*  -  *  . 


u-.u-.TT_i,  i-ir*  •  (-tr)~1~ 

'-^(w-£  -  ■-■•••' -(2-io) 


l/*4 1  =  —  (C/’+t/p) 


where 

A-- diag[<l-/9i)  Ifi'  A*  =  dia*[(l+^,)  i«] 

.4,-diag  (0,  4J),  0,  ^>o  (A) 

B  is  a  constant  to  be  determined.  (aU/ax)^*  and  (aU/ax)^ 
represent  the  first  order  backward  and  forward  difference 
quotient  at  point  j,  respectively.  Let  F  =  +  F^.  By 

rearranging  the  above  formula  it  becomes 


uy'=u:  +-£-[<(/;♦*- urri)+  (UV1  -£/■>] 


(2.11) 


because 


’  -(tc-X:  ■  — et  «"**«  -  »:*r» 


Let  d^=0^xj^>O,  then  Ag=diag  (d^).  Let 

n  _  c  .  c  ~  1 


D  =  S  A„S 

p 


(2.12) 


(2.13) 


then  to  a  second  order  accuracy,  equation  (2.11)  can  be  written 
as  the  following: 


(2.14) 


{I  +  D>  ^)dUV~,= 

(/  -f-  or*  &)x/ri —  *L(Fjn-Fr 

u:"=u;  ^-Lauf^-dur1) 

« 


i) 


n« » ••  ALai? 
U'"  A**' 


■  - 1 
i  ♦ » 


where  I  is  a  unit  matrix. 

In  the  following  let  us  investigate  the  problem  regarding 
the  selection  of  d^.  The  d^  which  satisfies  the  stability 
condition  is: 


d,^  max 


{[-r(^+£H£j.  •} 


(2.15) 


(1-1,  2, 


■n) 


Therefore,  we  may  choose 


{[A  (Uii+|l)^],0}  C216) 

This  is  MacCormack’s  method.  However,  just  as  we  pointed  out  in 
the  introduction,  this  choice  requires  operation  involving  two 
diagonal  matrices.  In  order  to  further  simplify  the  computation 
d^ ,  may  be  chosen  as : 


. d.-mtx  •  ©}  =  <*  (2.17) 

Here,  d  is  not  varying  with  1.  It  satisfies  the  largest  among  n 
parameters  in  the  stability  condition.  Substituting  equation 
(2.17)  into  (2.13)  we  get 
D  =  dl 

By  substituting  this  equation  into  (2.14)  we  get: 


(1  +  dj 71  |£)*/:F -  -  <*7n  -  FT77)  +  ^TTI  wvi 

t/:*,  =  c/j+-|-(3£/rT+«/TTT) 


(2.18) 


This  is  a  simple  implicit  calculation  scheme  consisting  of  an 


scalar  equations.  Furthermore,  each  scalar  equation  is 


independent.  They  can  be  solved  by  the  scalar  chasing  method. 


In  addition,  it  is  also  possible  to  change  equation  (2.14)  into  a 


characteristic  form.  In  that  case,  each  equation  will  have  its 


own  one-to-one  corresponding  d^  (1  =  1,2, . n).  By  using 

these  local  dj J s ,  we  also  can  obtain  an  implicit  formula  with  n 


scalar  equations.  Obviously,  that  method  is  also  very  simple. 


3.  Example 


The  sample  implicit  method  described  above  was  used  to 


calculate  the  laminar  separation  of  a  supersonic  flow  around  a 


two-dimensional  compressible  corner.  The  incoming  flow  condi¬ 


tions  are:  M  =  3.0,  ReB  =  1.68  x  10  ,  a  =  10°.  The  wall  is 


adiabatic.  The  initial  N-S  equation,  coordinate  transformation 


and  lattice  formation  are  identical  to  those  used  in  reference 


[2].  In  the  calculation  rigorous  implicit  condition  was  used  for 


the  body  surface  boundary.  Figures  1  and  2  show  the  distribution 


surface  pressure  and  that  of  frictional  drag,  respectively. 


Furthermore,  the  results  are  also  compared  with  those  obtained 


with  the  explicit  method  in  reference  [2]  and  those  with  im¬ 


plicit  formula  in  reference  [1].  Results  show  that  the  time 


for  each  step  is  less  than  that  of  MacCormack's  method.  But, 


*  *  -  \T  V  "A’  v  -.V -■W  V  •/  v  -/  v  %•  /  •/  •. 
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the  overall  efficiency  is  comparable.  Similar  to  the  situation 
in  reference  Cl],  after  the  computation  is  done  to  a  certain 
stage,  the  time  step  must  be  continuously  reduced  to  improve 
accuracy. 


o*  !•* , 

Figure  1.  Distribution  of  Wall  Surface  Pressure 
o  explicit  method 

-  maximum  stability  condition  method 
a  method  used  in  reference  [1] 


G.iW 


••ir , 

Figure  2.  Distribution  of  Wall  Surface  Frictional  Drag 


o  explicit  method 

-  maximum  stability  condition  method 
A  method  used  in  reference  [1] 

Comrade  Jie  Xinxing  also  participated  in  a  portion  of  this 

work . 
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A  SIMPLER  IMPLICIT  METHOD  SOLVING  N-S  EQUATION 
Zhang  Hanxin,  Ya  Zechu,  Lu  Linsheng 

<Tk«  Chian  AtrrffUMici  JlwMrcfc  and  Development  Caure) 

Abstract 

,  i  *  -  •  .  4  ■ 

A  simpler  implicit  method  based  oa  the  reference  Cl3  is  presented. 
The  process  of  iteration  is  .on)/  of  scalar  operation.  The  calculation  shows 
that  the  results  are  in  agreement  with  those  given  by  the  explicit  method, 
and  that  the  method  is  as  efficient  as  MacCormack's,  but  present  method 
is  simpler. 


Numerical  Calculation  of  Separation  Flow  Over  Severely  / 

Indented  Blunt  Body 
Gao  Shuchun  and  Gu  Gangmin 
(China  Aerodynamic  Research  and  Development  Center) 

During  re-entry,  a  vehicle  experiences  mostly  turbulent  flow 

in  the  process.  Due  to  high  temperature  chemical  attack  and 

particle  cloud  erosion,  the  head  is  indented.  Thus,  the 

calculation  of  separation  flow  over  an  indented  body  is  very 

significant.  Due  to  flow  separation  and  some  complicated  flow 

field  interference  effects,  the  structure  of  the  flow  field 

becomes  very  complex.  In  recent  years,  although  this  type  of 

f  1-31 

research  has  been  initiated  abroad  ,  yet  there  is  very  little 

progress.  In  particular,  calculated  results  of  turbulent  flows 

are  far  apart  from  experimental  data.  This  is  primarily  related 

to  factors  such  as  the  selection  of  the  turbulence  model,  the 

determination  of  the  position  of  the  turning  point,  and  the 

roughness  of  the  surface.  In  particular,  the  mechanism  of  the 

turbulence  is  still  not  very  clear.  Consequently,  there  is  some 

blindness  in  the  determination  of  models  and  the  establishment  of 

numerical  methods.  Many  people  doubt  the  validity  of  using  the 

hybrid  theory  to  study  the  separated  turbulent  flow.  Despite 

these  reasons,  people  are  still  using  the  simplest  turbulent  flow 

[1  2] 

model  to  simulate  a  flow  field  with  separation  present  *  .  We 

believe  that  this  is  the  simplest  treatment  before  many  compli¬ 
cated  phenomena  are  understood.  Moreover,  it  is  able  to  solve 
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some  practical  problems. 


In  this  work,  the  simplified  NS  equation  was  used  as  the 

starting  point  to  calculate  hypersonic  flow  field  over  a  severely 

indented  axi-symmetric  blunt  body  by  using  a  time  dependent 

r  i  ] 

implicit  spatial  factor  decomposition  method  and  by  adopting 

r  o  /i 

the  Cebeci-Smith  double  layer  vortex  viscosity  approximation1  ’ 
The  results  are  more  or  less  in  agreement  with  the  calculated  and 
experimental  results  obtained  abroad. 

The  example  used  in  this  work  is  Re  =6  x  10^  M  =9.0,  T  =200 

r./r.= 3.44,  •  •  08 

R,y^and  the  shape  of  the  object  is  the  same  as  that  in  reference 
[23. 

Figure  1  is  the  pressure  distribution  in  the  wall  surface, 
which  is  close  to  those  described  in  references  [1,23. 


Figure  1 

1.  this  work 
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V-V-\'V 

V  *1  ^  . 


Figure  2 


Figure  2  is  the  heat  flow  distribution  on  the  wall  surface 
which  is  close  to  the  result  reported  in  reference  Cl].  It  is 
also  close  to  the  experimental  value  on  a  small  wall.  However, 
it  is  quite  different  from  the  result  reported  in  reference  [2]. 
There  are  many  reasons  for  this  discrepancy.  For  instance,  in 
the  initial  equation  the  spacing  between  meshes,  the  numerical 
method,  and  the  specific  algebraic  model  can  lead  to  such 
differences.  We  believe  that  these  factors  are  not  essential. 

Why  is  there  such  a  large  discrepancy  in  the  heat  flow  when  the 
results  of  wall  pressure  are  basically  the  same?  This  is  because 
the  heat  flow  is  not  a  directly  calculated  quantity.  Instead  it 
is  a  quantity  derived  from  the  calculated  result.  Its  magnitude 
is  determined  by  the  temperature  gradient  and  thermal 
conductivity  near  the  wall  surface: 

This  paper  was  received  on  June  10,  1984. 


Figure  2 


Figure  2  is  the  heat  flow  distribution  on  the  wall  surface 
which  is  close  to  the  result  reported  in  reference  Cl],  It  is 
also  close  to  the  experimental  value  on  a  small  wall.  However, 
it  is  quite  different  from  the  result  reported  in  reference  [2], 
There  are  many  reasons  for  this  discrepancy.  For  instance,  in 
the  initial  equation  the  spacing  between  meshes,  the  numerical 
method,  and  the  specific  algebraic  model  can  lead  to  such 
differences.  We  believe  that  these  factors  are  not  essential. 
Why  is  there  such  a  large  discrepancy  in  the  heat  flow  when  the 
results  of  wall  pressure  are  basically  the  same?  This  is  becaus 
the  heat  flow  is  not  a  directly  calculated  quantity.  Instead  it 
is  a  quantity  derived  from  the  calculated  result.  Its  magnitude 
is  determined  by  the  temperature  gradient  and  thermal 
conductivity  near  the  wall  surface: 


134 


Figure  3  compares  the  frictional  drag  distribution  in  a 
laminar  flow  to  that  in  a  turbulent  flow.  The  frictional  drag  in 
the  turbulent  flow  was  calculated  using  a  formula  similar  to 
equation  (3). 

Figure  4  compares  the  shock  wave  position  in  a  laminar  flow 
to  that  in  a  turbulent  flow,  which  agrees  with  the  conclusion  in 
reference  [23. 

Figures  5-7  show  the  iso-M,  iso-density  and  isobaric  contours 
of  the  turbulent  flow,  respectively.  The  structure  of  the  flow 
field  is  described  in  detail. 

From  Figure  5  we  can  see  that  the  shock  wave  becomes  steeper 
after  the-  separated  flow  is  re-attached.  Hence,  there  is  a 
second  subsonic  region. 


Figure  3 

1.  turbulent  flow 

2.  laminar  flow 


Figure  4  Comparison  of  Laminar  and  Turbulent  Shock  Wave 

1.  turbulent  flow 
2..  laminar  flow 


Figure  5  Iso-Mach  Number  Lines  for  Turbulent  Flow 

Figure  8  shows  the  flow  pattern  in  the  separated  region. 
There  is  a  "Secondary  Separation  Bubble",  which  is  in  agreement 
with  the  conclusion  in  reference  [5].  This  effect  was  not  found 
in  calculating  a  turbulent  flow. 


After  approximately  3000  iterations,  the  flow  field  of  the 
turbulent  flow  could  basically  be  created.  We  were  inspired  by 


the  work  in  reference  [2]  to  observe  this  effect.  Our  conclusion 
also  agrees  with  theirs.  The  authors  of  that  paper  believed  that 
it  was  caused  by  the  interference  of  the  separation  bubble  and 
the  subsonic  band  of  the  blunt  body.  They  also  pointed  out  that 
a  similar  effect  could  be  found  in  calculating  a  laminar  flow. 
However,  we  did  not  observe  this  oscillating  phenomenon  in  our 
laminar  flow  calculation. 


Figure  6  Iso-density  Contours  of  Turbulent  Flow 
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Figure  7  Isobaric  Lines  of  Turbulent  Flow 


Figure  8  Velocity  Trend  in  Separated  Laminar  Flow 

With  a  time-dependent  technique  to  solve  a  separated  flow 
field,  because  the  variation  rate  of  a  physical  quantity  with 
time  oscillates  periodically,  the  flow  field  does  not  converge 
consistently.  However,  in  a  certain  sense  of  averaging,  the  flow 
field  is  convergent.  Therefore,  we  should  establish  specific 
criteria  for  an  average  convergence.  Of  course,  such  criteria 
may  vary.  The  NS  equation  for  a  laminar  flow  is  quasi-linear 
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and  the  turbulent  flow  equation,  even  in  the  simplest  algebraic 
model,  is  non-linear.  There  are  difficulties  in  the  solution 
finding  process  which  are  not  encountered  in  the  laminar  flow 
calculation.  Therefore,  we  cannot  copy  a  laminar  flow  model  in 
the  calculation  of  a  turbulent  flow. 
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NUMERICAL  CALCULATION  OF  SEPARATION  FLOW 

OVER  SEVERELY  INDENTED  BLUNT  BODY 

*  - 

Gao  Shuchun,  Gu  Gangmm.  ...  . 

(China  Aerodynamic  Reitareh  tt  Development  Centre ) 

\  '  \ 

,  Abstract 

A  time  dependent  space  factored  implicit  numerical  procedure  solving 
simplified  Navier-Stokes  equations  is  used  to  calculate  axi-sy mmetric 
hypersonic.  Tiseous  flow  over  a  severely  indented  blunt  body.  A  simple 
algebraic  turbulent  modal  is  adopted.  The  results  are  compared  with  ano¬ 
ther  numerical  solution  in  ref  1.  A  good  agreement  is  found  for  the  wal 
presure  and  heating  distribution,  although  there  are  some  discrepancy  m 
detail  for  flow  structure.  A  secondary  separation  bubble  is  seen  clearly  m 
laminar  numerical  results  and  high  frequency  oscillation  phenomena  are 
found  to  be  the  same  as  that  in  the  turbulent  calculation  of  ref  2. 


An  Implicit  Technique  for  Calculating  Separated  Base  Flow 

Dong  Changquan 

(Beijing  Institute  of  Information  and  Control) 

1.  Introduction 

In  solving  the  Navier-Stokes  equation  with  an  explicit 
technique,  its  stability  is  often  limited  by  the  time  steps. 

When  it  is  applied  to  a  high  Reynolds  number  viscous  flow,  we 
have  to  pay  a  high  price  even  in  a  two-dimensional  case.  One  of 
the  major  problems  in  any  implicit  technique  is  the  treatment 
of  the  non-linear  convection  and  pressure  terms.  Earlier  com¬ 
putation  'shows  that  an  alternating  direction  implicit  scheme  in 
conjunction  with  the  quasi-linearization  of  non-linear  terms 
could  be  used  to  calculate  the  separated  base  flow  field.  This 
type  of  scheme  is  obviously  able  to  accelerate  the  stabilization 
process,  as  compared  to  explicit  method,  in  order  to  save  compu¬ 
ter  time. 

The  scheme  used  in  the  work  is  based  on  the  SCM  technique 
developed  by  Steger  and  Warming  for  solving  Euler  equations  and 
MacCormack’s  1980  algorithm.  (The  advantages  of  these  schemes 
is  to  avoid  solving  three  diagonal  or  three  diagonal  block 
equations  and  to  directly  create  an  iteration  scheme  to  simplify 
the  mathematics  in  each  iteration.)  A  matrix  similarity  trans¬ 
formation  is  used  to  solve  the  N-S  equation  by  the  SCM  technique 

r  6 1 

The  viscosity  terms  are  treated  by  CayflbeB's  treatment  for 
diffusion  terms  as  well  as  by  a  technique  similar  to  that  used 


by  Allen-Cheng-Bpau*oBcKa5U  This  algorithm  can  eventually  derive 
an  upper  or  lower  triangular  matrix.  The  direct  iteration  method 
can  save  the  amount  of  computation  in  each  iteration  as  compared 
to  a  tridiagonal  chasing  method. 

Finally,  it  should  be  pointed  out  that  in  order  to  meet  the 
stability  requirement  it  is  very  important  to  add  an  appropriate 
dissipation  factor  in  the  process.  In  addition,  the  boundary 
condition  must  be  approached  in  a  way  consistent  with  that  in  the 
SCM  technique.  The  algorithm  has  an  accuracy  to  the  second 
order . 

2.  Basic  Equations 

-  The  following  is  the  non-dimensional  Navier-Stokes  equation: 

^dF.  dc»_  d(y,  +  y,)  dtiy.+jy,)  (1) 

dt  dx  dy  dx  dy  ~  di  dy 

The  reason  to  decompose  the  viscosity  terms  and  5^  into  +V2 
and  V?! +W2  is  to  facilitate  the  treatment  of  mixed  derivatives. 

It  also  facilitates  the  use  of  CaynbeB's  or  Allen-Cheng-Bpau*- 
oBcKafl's  viscosity  treatment.  A  detailed  expression  is  shown  in 
reference  C4],  The  region  of  computation  and  the  boundary 
condition  are  defined  using  the  same  method  as  that  used  in 
reference  [4] . 

3.  Computation  Technique 

Almost  at  the  same  time  as  MacCormack  introduced  his 

r  3 1 

algorithm  in  1980,  Steger-Warming  presented  the  SCM  method 
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for  solving  Euler  equations.  Furthermore,  based  on  whether  the 
eigen  value  of  the  Jacobi  matrix  A  or  B  is  positive  or  negative, 
the  derivatives  are  approximated  by  off-center  differences.  The 
coefficient  matrix  of  the  derived  linear  algebraic  equation  is  an 
upper  or  a  lower  triangular  matrix  to  facilitate  direct 
computation.  Consequently,  the  tridiagonal  chasing  method 
commonly  required  in  the  Beam-Warming  factoring  method  can  he 
avoided  to  reduce  the  load  of  computation  in  the  stabilization 
process.  We  employed  this  technique  to  solve  the  Navier-Stokes 
equation  with  the  objective  of  making  the  computation  more 
economical . 

_  According  to  the  second  order  three  layer  format  introduced 
f  31 

by  Beam-Warming  : 


Ay'1 


0A# 


±  <^-,+7Tr^r+TfsA‘?- 


l  +  £  dt 

From  equation  (1)  we  get 

&0'+iTe{^A+-kB)  A°--nk{i£r+&y+*f,s 


where 


B- 


&u  *u  &u.  &U  dU, 


(2) 


Because  F  and  G  are  homogeneous  functions  of  U,  it  is  possible  to 
introduce  the  following  similarity  transformation: 


/80 


Therefore : 


where 


i '  +  TFT  (£' ^ +  £ *-■ S* +£  B')} 

F'  =  A*0  F'-A-O  G*  =  B+u  G~  =  B~U 


F  (At»  A«)  — - — 

2y 


2  (y  —  1)  Aj + Aj  +  A< 

2(y-l)A,ti+Al(ii+c)+A»(ii-c) 

2(y-l)A,y+A,ti+y4w 

(y-DA.W+yo+A  C(*+c)*+o*]  + 
2 


+  AC(#_c).  +  ^3  +  WrI 

>  V.  '  a  •„  ••  ■ 

i^la  (3— y) (A,+A«)c* 

2(r~i> 


(  2(y-l>  Aj+  A>  +  A  • 

„  p  1  2<y—  1)  £  ,u  +  £  ,u+  x  ,u 

G  <  A  ■  >  A  >  A J  2y  2<y- 1)  ,y  +  £  ,(v  +  c)  +  £  ,(v-c) 

i  (y-i>5,  (u‘  +  y‘)  +Ai  C«‘+ (y+c)*]  + 

w  <3-y)(I,+  3«)c» 

2(y- 1) 

F±  =  F  (At,  Af,  AJ) 

C  =G  <  A  i»  A  *»  A») 

The  splitting  of  the  eigen  value  is  not  unique.  There  are  other 
ways  to  divide  it  into  positive  and  negative  parts.  In  order  to 
overcome  the  critical  state  of  \+  and  A-,  appropriate  small 
quantities  Ex  and  e2  are  introduced.  Hence 
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F*  x  F  (A*  +  eit  Aj'-r c,,Al  +  i,),  F  =F(A,  -e,,  A,  e,,  At  —  et) 

G*  =(?(  4  I +  4i+<»»  A  I '♦'*«)»  ^  =  G(  ^  i  —  Cn  A»— e:>  A*- 


In  the  following,  we  will  discuss  the  modification  of  the  left 
hand  side  terms  by  RHS.  Let  us  approach  second  order  terms  by 
CaynbeB's  method''^: 


dx‘  A*1 
d‘F'*>  l  , 


■♦i  i  ,  _ 

^ - =  (  R‘  •  J  _  ©*  *  's  _  1  /  p*  ♦  i  n«  *  i . 

dx1  AxJ  *  a,i 


The  computation  format  after  the  modification  is: 

[  1  V.(^-  +  >4*)]  [  I  +  (1*BV  *■«*->]  4(?* 

*  WF*  +  6iG*+VF'-  +  3’G-)  +  RtfS 

[ 1  +  7T+1W  v.<fl-+fi*)][i  +  Y'lf'-j,  4.or > ] 4y ^ . 4 1? ' 

t/^rri=  u'  +  &u~l 

L  £  £  L 

6°  6“.,  6r  and  6°  are  second  order  single  side  difference 

x  x  y  y 

operat°rS  •  [ 1  + ( 1  + 

-^T<VF*+*lG*+if.F’  +  *lG')  +  RHS* 


(1+1)  Ay 


a, (*-+*•>]  [I’t,Ti+f)srA,<ir)]  *0'"-*$' 


(.u'+0'*'+ A(7,4*> 

Z 

It  should  be  pointed  out  that  a  filtering  technique  consistent 
with  the  algorithm  is  required  for  regions  with  intense 
interruption  and  high  gradient  such  as  shock  waves.  In  order  to 
satisfy  the  stability  requirement,  it  is  very  important  to 
include  some  fourth  order  dissipation  factors  in  the  process  such 
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a>.  Ajc*  a4Lf  ,  A y‘  d'Q 

8  (l  +  £)  dx‘  a ’  (1  +  |)  ay4 

u>  and  to  are  the  dissipation  coefficients.  The  entire  scheme  is 
x  y 

[2] 

accurate  to  the  second  order.  Just  as  Beam-Warming-Yee 
pointed  out,  the  appropriate  approximation  of  the  implicit 
boundary  condition  is  very  important  in  an  implicit  scheme.  It 
is  able  to  createan  unconditionally  stable  algorithm  consistent 
with  the  format  of  the  internal  points.  (Obviously,  the  boundary 
must  exist  for  (At/Ax)  and  (At/Ay)  when  At,  Ax  and  Ay-  0.)  Due 
to  page  limitation,  this  cannot  be  discussed  here. 

4.  Computational  Results 

In  this  work,  an  attempt  was  made  to  derive  a  direct,  non¬ 
iterative  implicit  scheme  to  overcome  the  deficiency  of  the 
traditional  upwind  format.  The  resolution  of  the  flow  field  was 
improved  because  the  off-center  approximation  was  based  on  the 
sign  of  the  eigen  value  of  the  coefficient  matrices  A  and  B. 

We  calculated  the  results  in  several  cases  corresponding  to 
and  Reoo=4000,  1.5  x  10^  on  a  Univac  1100  computer. 

Compared  to  the  method  used  in  reference  [4],  each  iteration 
could  save  about  a  half  of  the  computer  time.  However,  it 
requires  more  iterations  to  converge.  The  entire  calculation 
saves  approximately  20%  computer  time  as  compared  to  the  process 
used  in  reference  [4],  The  results  obtained  with  both  methods 
are  basically  consistent.  The  velocity  distribution,  density 
distribution  and  temperature  distribution  on  the  x=1.5  cross- 
section  are  shown  below,  in  comparison  to  the  results  reported 
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in  reference  [4], 


Figure  1 


Figure  2 


Distribution  of  Velocities  u  and  v  on  the  x=1.5 
Cross-section 

1.  this  work 


Density  Distribution  p  and  Temperature  Distribution 
T  on  the  x=1.5  Cross-section 


1 .  this  work 
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AN  IMPLICIT  TECHNIQUE  FOR  COMPUTATION 
OF  BASE  FLOWFIELD 

Dons  Changquan 

(Beijing  Institute  of  Injormatlan  ft  Control) 

Abstract 

A  aew  noniterative  implicit  scheme  for  computation  of  the  base  flow- 
B4Ma  *ro  presented.  The  new  implicit  scheme  for  solution  of  Navier-Stokes 
eqs.  based  on  the  SCM  technique  is  developed  by  Steger  &  Warming  (I) 
for  solution  of  Euler  eqs.  and  some  properties  of  the  algorithm  analogous 
to  the  MacCormack  "algorithm  (5). 

Unlike  usual  implicit  schemes,  the  present  algorithm  does  not  require 
the  solution  of  a  block  tridiagonal  systems,  because  it  leads  to  block  bidia¬ 
gonal  systems.  Numerical  results  indicate  rather  improvement  in  accuracy 
for  circulation-region  than  (4).  Futhermore  the  computer  time  required  to 
obtain  a  converged  solution  has  been  reduced. 
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Computation  of  Compressible  Turbulent  Separated  Boundary  Layer 

Liu  Songling 

(Northwestern  Polytechnical  University) 

I. 

It  has  always  been  an  important  problem  in  fluid  dynamics  to 

predict  a  separated  flow.  Although  the  N-S  equation  can  describe 

a  separated  flow,  however,  it  requires  a  great  deal  of  computer 

time  to  solve  this  equation  using  available  techniques.  Usually 

a  method  to  solve  the  boundary  layer  equation  fails  due  to  the 

presence  of  singular  points  near  the  separation  point.  In  recent 

years,  inverse  methods  to  solve  the  boundary  layer  equation, 

i.e.,  to  determine  the  boundary  velocity  u  outside  the  boundary 

layer  and  other  parameters  with  given  boundary  layer  displacement 

thickness  6*  or  friction  coefficient  c^,  have  been  developed. 

Since  both  displacement  thickness  and  friction  coefficient  are 

unknown,  it  is  necessary  to  first  estimate  these  parameters  in 

using  the  inverse  method.  The  calculated  u  value  is  then 

e 

compared  to  that  of  an  inviscous  flow  to  modify  6*  or  c^  for 
recalculation  until  the  velocity  distribution  obtained  agrees 
with  that  of  inviscous  flow. 

If  a  permeation  equation  is  used  as  an  auxiliary  equation, 
then  the  incompressible  flow  equation  formed  by  this  auxiliary 
equation  and  the  momentum  integral  equation  is 


(1) 


dH 

dx 
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0  is  the  boundary  layer  momentum  thickness  and  H  is  the  shape 
factor.  H,=(6-6  )/9,6  is  the  boundary  layer  thickness,  and  Cg  is 
the  permeation  coefficient.  When  a  forward  mode  is  used,  (1/u^) 
(duc/dx)  is  known,  c^  and  Cg  can  be  calculated  empirically,  and 
dH/dH,  can  be  determined  from  ^=^00  which  has  already  been 
established.  Equations  (1)  and  (2)  can  be  written  as: 
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As  H  increases,  H,  decreases  and  dH,,  /dH  approaches  zero.  There 
is  a  singular  point  in  equation  (3).  If  6*  and  uc  are  specified 
as  the  knowns ,  then  equation  (3)  becomes: 
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Because  dH, /dH  <0,  the  determinant  of  the  coefficients  on  the 
left  side  is  always  not  equal  to  zero  to  avoid  any  singular 
point. 

[  1  ] 

The  Lag-Entrainment  method  developed  by  Green  et  al  is 
widely  used  in  the  computation  in  aerodynamics.  This  material 
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was  used  in  reference  [8]  to  calculate  the  unseparated  turbulent 
boundary  layer  on  the  surface  of  a  compressor  blade  with 
satisfactory  results.  In  reference  [2],  the  shape  factor 
relation  Hx  =HX  (H)  was  modified  to  be  used  in  an  inverse  mode 
calculation  in  an  example  of  separation  caused  by  a  shock  wave. 

In  this  work,  this  method  is  used  to  calculate  the  separated 
compressible  flow  in  the  trailing  edge  by  using  the  measured 
displacement  thickness  to  determine  the  flow  velocity  ue, 
momentum  thickness  e  and  friction  coefficient  c^  on  the  boundary 
layer . 

From  the  momentum  integral  equation  and  permeation  equation 
we  can  derive  that 


0  du,  _  1  /  d<3*  p 
v,  dx  Ft\dx 


+  (-1  +  0 

2 


Hlc,\  djj  ' 
2  JdH, 


dfj 


F,=  —  H (.H  +2  —  M l)  +  (1  +  0 ,2rA/J)  (//  +  !)//, 


(5) 


In  addition,  H=6*/0  and  H=(H+1 ) / ( 1+0 . 2rM* )-1 ,  where  r  is  the 
restoration  coefficient.  The  relation  between  and  ff  is  the 
one  recommended  in  reference  [2]  .  Equation  (5)  and  reference  [1] 
give  a  closed  set  of  three  differential  equations.  After  the 
initial  values  re  specified,  ug,  H,  0  and  c^  can  be  solved. 

It  is  demonstrated  in  practice  that,  when  the  initial  point 
is  in  an  area  with  a  larger  shape  factor,  the  value  of  c^ 
obtained  using  the  method  recommended  in  reference  [1]  may  be 
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too  large,  leading  to  small  ug  and  H  along  the  process  (Example 
2).  After  changing  to  the  method  used  in  reference  [3]  to 
determine  the  initial  value  of  Cg  and  using  the  same  differential 
equation  as  the  one  in  reference  [1]  to  calculate  Cg  along  the 
process,  the  result  is  improved  significantly. 


The  method  described  above  was  used  in  the  following 
examples : 

(1)  Boundary  Layer  on  the  Upper  Surface  of  a  Supercritical 
Airfoil . 
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Figure  1 
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Figure  2 

1 .  this  work 

2.  experimental 

Figure  1  shows  the  measured  velocity  distribution  and 
boundary  layer  on -an  airfoil  as  reported  in  reference  [4],  The 
incident  flow  M  number  is  0.75  and  Re  number  is  2  x  10^.  It  is 
separated  near  the  trailing  edge.  The  figure  shows  calculated 
results  obtained  by  using  forward  schemes  introduced  by  Cebeci- 
Smith  and  Nash-MacDonald .  In  the  range  x/c<  0.85,  a  forward  mode  / 86 
was  used.  There  is  a  shock  wave  at  x/c  =  0.2  to  make  the 
boundary  layer  thickness  grow.  But,  it  is  not  separated.  At  x/c 
=0.85,  it  was  switched  to  an  inverse  mode.  The  displacement 
thickness,  momentum  thickness  and  friction  coefficient  are  in 
agreement  with  the  experimental  data.  The  calculated  result  at 


the  trailing  edge  is  much  better  than  those  obtained  with  both 
forward  modes.  The  calculated  M  number  at  the  external  boundary 
of  the  trailing  edge  boundary  layer  is  less  than  47«  different 
from  the  experimental  value. 

(2)  Boundary  Layer  on  the  Upper  Surface  of  a  Circular  Arc 
Airfoil  in  Reference  [5] 

The  thickness  of  the  airfoil  is  107o,  the  chord  length  is  12 
inches,  the  Re  number  is  3.6  x  10^,  and  the  incident  M  number  is 
0.7339.  There  is  a  shock  wave  at  0.62  chord  length.  The  M 
number  in  front  of  the  shock  wave  is  1.32,  causing  a  turbulent 
boundary  layer  separation.  At  approximately  1.05  chord  length, 
the  boundary  layer  is  readhered.  Figure  2  shows  the  calculated 
results,  as  well  as  the  experimental  results  given  in  reference 
[53.  The  calculated  M  number  and  shape  factor  are  in  good 
agreement  with  the  experimental  data.  The  calculated  friction 
coefficient  is  slightly  low  in  the  separated  region.  This  method 
can  be  used  to  predict  the  separation  and  re-adhesion  of  the 
boundary  layer.  The  dotted  line  in  the  figure  indicates  the 
result  obtained  using  the  initial  value  of  c£  as  calculated  by 
the  method  adopted  in  reference  [13.  After  changing  to  the 
method  used  in  reference  [33  to  calculate  the  initial  value  of 
Cg,  the  solid  line  in  the  figure  is  obtained. 

(3)  Boundary  Layer  on  the  Upper  Surface  of  a  Transonic 
Compressor  Blade 

The  experimental  data  obtained  on  the  upper  surface  boundary 
layer  of  the  transonic  compressor  turbine  blade  is  reported  in 


reference  161.  The  upstream  M  number  is  0.85,  chord  length  is 

0.08  M  and  Re  number  is  5  x  10^.  The  calculated  results  are 

shown  in  Figure  3.  As  compared  to  the  experimental  data,  the 

maximum  deviation  of  M  is  less  than  5%.  The  momentum  thickness 

e 

is  on  the  low  side.  After  (x/c)>  0.75,  c^  is  negative.  Thus,  we 
believe  that  the  turbulent  boundary  layer  is  already  separated  in 
this  region. 

(4)  Boundary  Layer  Reported  in  Reference  [7] 

The  experimental  specimen  is  a  circular  arc  with  a  thickness 

of  18%.  The  incident  M  number  is  0.7  and  the  Re  number  is  4x10^. 

Figure  4  shows  the  results  obtained  by  using  the  method  described 

in  this  paper,  together  with  the  solution  of  the  N-S  equation 

given  in  reference  [7].  The  M  number  at  the  trailing  edge 
calculated  in  this  method  is  87*  lower  than  the  experimental 
value.  The  momentum  thickness  is  slightly  larger  than  that 
calculated  based  on  the  N-S  equation.  Both  methods  agree  with 
the  experimental  data  to  the  same  extent.  The  measured 
separation  point  is  located  at  x  =  -0.02M.  The  calculated 
separation  point  using  either  method  is  located  downstream  from 
the  measured  point.  Based  on  this  example,  we  find  that  the 
inverse  mode  technique -is  not  very  different  from  solving  the  N-S 
equation. 


uation  3.  experimental 


56 


III. 
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From  the  above  computation  we  arrive  at  the  following 
conclusions : 

(1)  It  is  possible  to  avoid  the  singular  point  near  the 
separation  point  when  using  an  inverse  method  to  solve  the 
boundary  layer  equation.  It  can  be  used  to  calculate  a  separated 
boundary  layer  on  the  trailing  edge  when  the  intensity  is  not  too 
high.  The  method  described  can  be  used  to  calculate  turbulent 
separated  boundary  layers  with  satisfactory  results. 

(2)  Examples  show  that  the  result  of  the  inverse  method  near 
the  separation  point  is  not  too  different  from  that  obtained  by 
solving  the  N-S  equation. 
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THE  COMPUTATIONS  OF  COMPRESSIBLE  TURBULENT 
SEPARATED  BOUNDARY  LAYERS 

t 

Liu  Songling 

(The  North-western  Polytechnical  University) 

Abstract 

The  lag-entrainment  method  is  used  in  an  inverse  mode  to  predict  two 
dimensional  compressible  turbulent  separated  boundary  layers  in  the  trai¬ 
ling  edge  region.  There  are  reasonable  agreements  of  the  predicted  values 
#irh  those  several  measured  from  separated  turbulent  boundary  layer*. 
As  example  shows  that'  the  differences  between  the  results  obtained  by 
present  method  and  the  N-S  equation  solutions  are  not  significant. 
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An  Inverse  Boundary  Layer  Method  for  Separated  Flow 

Qin  Ning 

(Nanjing  Aeronautical  Institute) 
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The  inverse  boundary  layer  method  is  a  numerical  method  for 
solving  boundary  layer  flows.  In  this  method,  the  displacement 
thickness  6*  or  the  wall  shear  stress  is  the  prescribed 
variable  and  the  velocity  distribution  ug  or  the  pressure 
distribution  of  the  boundary  edge  flow  is  calculated  with  the 
boundary  layer  solution.  This  can  avoid  the  singularity  at  the 
separation  point  in  the  finite-difference  calculation  and  provide 
a  meaningful  method  in  solving  separated  flows.  This  method 
requires  far  less  calculations  and  storage  space  compared  to  the 
direct  solution  of  the  N-S  equation.  In  practice,  it  is  a  good 
method  of  approximation  for  solving  bubble  separation  in  the  boundary 
layer. 

We  have  based  on  the  inverse  boundary  layer  method  for 

separated  flows  given  by  Cebeci  et  al  to  propose  an  inverse 

T  21 

calculation  method  for  the  incompressible  boundary  flows  .  In 
this  method  the  displacement  thickness  6*  is  prescribed  and  the 
profile  of  the  boundary  layer  edge  velocity  is  calculated  with 
the  boundary  layer  solution.  Falkner-Skan  coordinates  are  used 
both  in  direct  and  inverse  calculation  regions  in  order  to  avoid 
the  matching  difficulty  between  the  two  regions  as  that  occurred 
in  reference  [2]  which  have  used  Falkner-Skan  coordinates  for  the 
direct  calculation  region  and  the  physical  coordinates  for  the 
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inverse  calculation.  Inadequate  treatment  of  the  solution  at  the 
transition  point  may  cause  jump  or  divergence  of  the  solution. 
With  the  transformation  coordinates,  the  variation  of  the 
solution  at  each  point  is  smaller  than  that  with  the  physical 
coordinates  and  the  convergency  is  faster  when  the  values  of 
previous  solution  are  used  as  the  initial  values  for  the  next 
calculation  in  Newton  iteration.  Additionally,  the  use  of  only 
one  transformation  coordinate  allows  the  use  of  one  general 
calculation  process  from  the  starting  point  x=0  of  the  direct 
calculation  region  to  the  inverse  calculation  region  containing 
separations . 

I.  Basic  Equations  and  Transformation 
The  two-dimensional  steady  state  boundary  layer  equations 
for  incompressible  laminar  flows  are  as  follows: 


(la) 


ox  oy 


du  du  .  du,  ,  c^u 
U  dx  V  dy  '  dx  dy  * 


(1b) 


V- o, 


if  — =  0  s 


y=v,f 


u  =  u. 


(2) 


In  direct  calculation,  ug  is  given. 


In  inverse  calculation,  ug  is  unknown  and  6*  is  given  and  the 


additional  boundary  condition  is  obtained  from  the  definition: 


Using  Falkner-Skan  transformation  for  the  above  equations 


and  boundary  conditions, 
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In  direct  calculation,  ufi  is  given. 

In  inverse  calculation,  ug  is  unknown,  and  the  boundary 
conditions  are 


y  vx 

where  f  is  the  dimensionless  stream  function,  <p=(uevx)% 
f'=u/u  ,  m=x/u  du  /dx,  the  superscript  is  a/an. 

6  v  6 


(7) 
f  (x,  n) , 


II.  Numerical  Method 
1.  The  Finite-Difference  Scheme 

The  Keller  box  scheme  is  employed.  It  has  the  second  order 
accuracy  with  unconditional  stability.  Multiplying  equation  (5) 
by  ue,  let  f=F  and  ug=W,  the  differential  equation  (5)  is  reduced 
to  the  first  order  differential  equations 

F'=U  (8a) 

U'=V  (8b) 

wv'  +  L  a  u-w>-.w(if£-i'4£-)  (8c) 


In  the  inverse  method, 

v=t,..  a* 


(10) 


where 

5.—^-.  X-j-,  i'-dWW^TlL 

With  X=X„  ,  set  up  an  equation  at  the  point  (X0  , n j ^  » 

with  X=X.(i>0),  set  up  an  equations  at  point  (X.jn^^j/a) 

of  the  two  previous  equations,  and  set  up  an  equation  at  point (X.-j. 
of  the  third  formula.  The  partial  derivatives  are  approximated 

using  first  order  central-difference  operators.  The  following 

finite-difference  equations  are  obtained. 


*7-» 

A;!,  — L(K,-K>.J)  =  o 
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and  the  boundary  conditions  are 
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F;=t/;=o, 


F'r 


(j*)1 

~Xin 


(12) 


The  solutions  of  3J+1  unknowns  (F j ,  U j ,  V i ,  W1)  can  be  obtained 
from  the  group  of  nonlinear  algebraic  equations  formed  from  the 
above  equations  and  the  boundary  conditions. 

2.  Linearization  of  the  Equations  and  Obtaining  the 


Solutions 

Using  Newton's  iteration  method  let 

F;'-F;+dFlt  U\=U\+SU,,  V\-V\+&V„ 
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and  substituting  them  in  equations  (11)  and  (12),  the  following 

2 

linearized  equations  are  obtained  with  the  omission  of  6  and 
other  small  quantities. 

h~,lt  (dF,  —  3F , - » )  — (S(J ,+dU =  (r,),  ,  ,  .  . 

2  (13a) 

h 7  i ,  (dU ,  -  dU , . , )  -  i-  (<3K ,  +  df' , . , )  =  (r , ) , . , 

(13b) 

(S.)  .dF,.,  -  (S,),dF,  4-  (5,)  +  (5.) ,«/,  -  (5S)  .ir,., 

^(5.)IaK)^(S,),6Jr=(r;), 

;  =  2,  3,  y  (13c) 

6F,=  (r,),  =  0,  dUi=(r.),  =  0 

<^F'+2()^F)T7^<J,K=(r,>,  (14) 

The  coefficients  form  a  4x4  diagonal  matrix  and  the  solution  of 
the  above  linear  equations  can  be  obtained  directly  using  the 
general  forward  method. 

The  converged  solution  of  a  fixed  point  Newton  iteration  can 
be  used  as  the  initial  value  to  obtain  the  solution  for  the  next 
point. 

3.  Treatment  of  the  Backflow  Region 

The  inverse  method  can  avoid  the  singularity  at  the 

separation  point.  However,  exponential  divergency  may  still 

occur  if  the  backflow  region  is  not  treated.  A  theoretical 

analysis  has  been  carried  out  in  reference  [6]  based  on  the 

suitability  of  the  initial  conditions  of  the  parabolic  equations. 

w«*  have  used  the  approximation  method  given  at  Reyhner  and 

[3] 

M  i*ge-[.otz  to  treat  this  problem  .  That  is 


f 'af '/9x=0 


at  u<0,  f,2=0, 

4.  The  Grids  of  Calculation  and  the  Growth  of  the  Boundary 

Layer 

With  transformation  coordinates,  the  thickness  of  the 

boundary  layer  of  the  attached  laminar  flow  is  basically 

constant.  Unevenly  distributed  grids  are  employed  along  the 

normal  direction.  The  grids  are  denser  near  the  wall.  Uneven 

grids  are  also  used  along  the  longitudinal  direction.  The  grids 

employed  in  the  direct  method  of  calculation  for  the  attached 

boundary  layer  flow  are  rather  sparse,  while  those  employed  in 

the  inverse  method  of  calculation  for  the  region  of  separation 

flow  are  denser.  Since  the  thickness  of  the  boundary  layer 

increases  due  to  separation,  the  grids  along  the  normal  direction  in 
the  separation  region  should  be  increased. 

III.  Example 

1.  For  testing  the  accuracy  of  the  numerical  method,  we 

*  t 

have  calculated  the  laminar  separation  flow  with  Re  T=10^,  and 

00  ,  L  ’ 

two  types  of  6*  distribution^ . 

The  calculation  started  from  X-1.0  and  the  Blasius  solution 
was  obtained.  Inverse  calculation  was  carried  out  starting  from 
the  second  point  with  given  distribution  of  6*,  aX=0.02.  There 
were  45  points  along  the  X  direction  and  41  points  along  the  X 
direction,  which  increased  to  61  points  after  the  separation 
point. 

The  results  of  the  numerical  calculation  are  shown  in  /91 

Figures  1  and  2  and  compared  with  the  results  given  by  Carter1-4-1 
using  global  iteration  calculations. 
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2.  Calculation  Based  on  Horton's  Ntitn  Experimental  Data 

Given  in  Reference  [5],  Direct  method  of  calculation  was  carried 

out  for  the  region  from  x=0  to  x=x^v  using  given  experimental 

values  u  .  Inverse  calculation  was  carried  out  for  the  region 
e 

starting  from  x=Xj^y  using  experimental  values  of  6*.  With  u# 
=9.36M/sec,  the  calculated  results  are  compared  with  the 
experimental  values  in  Figures  3  and  4. 


Figure  3. 


1 .  our  calculation 

2.  experimental  result  C 5 3 


Figure  4. 

1.  experimental  result  C 5 ] 

2.  -our  calculation 

3.  separation  point  in  experiment 

4.  separation  point  in  calculation 

The  calculation  was  carried  our  using  FACOM  230  computer. 

The  CPU  time  is  about  50  sec  per  calculation. 

IV.  Discussion 

By  comparison  of  the  results  of  calculation,  we  have 
concluded  as  follows: 

1.  The  calculation  results  using  the  forward  approximation 
method  of  Reyhner  and  Flugge-Lotz  are  consistent  with  the  results 
obtained  by  global  iteration  calculations.  However,  the  former 
method  used  less  calculations  and  storage  space  and  is  an 
effective  method  for  calculating  the  separation  of  bubbles. 

2.  .  Using  Falkner-Skan  coordinates  in  both  direct  and 
inverse  calculation  regions  can  provide  smooth  transition 
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between  the  two  regions  and  allows  the  use  of  the  same  numerical 
procedure  for  calculation.  No  jump  of  the  solution  will  occur  at 
the  transition  point. 

The  author  wishes  to  express  his  appreciation  to  Prof.  Cao 
Qi-peng  for  the  assistance  and  consultation  in  calculation  and 
the  preparation  of  the  manuscript. 
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Abstract 

This  paper  presents  an  numerical  method  for  solving  2-D,  steady,  in¬ 
compressible  laminer  boundary  layer  flows  with  separation  and  reattach- 
meat.  In  the  method,  the  displacement  thickness  is  prescribed  and  the  bo¬ 
undary  layer  edge  velocity  is  calculated  with  the  boundary  layer  solution 
so  that  the  singularity  at  separation  point  is  avoided,  in  the  reverse  flow 
region,  FLARE  approximation  is  used.  The  finite-difference  scheme  used 
here  is  2nd  order  accuracy  and  unconditionally  stable  implicit  Box  meth¬ 
od.  Falkner-Skan  coordinates  are  used  both  in  direct  and  inverse  calcula¬ 
tion  regions  in  order  to  avoid  the  matching  difficulty  between  the  two 
regions.  The  results  are  compared  with  global  iteration  calculations  and 
experimental  data. 


Wall  Lift  Interference  Corrections  in  Ground  Effect  Testing 

Li  Jingbai  Qi  Mengbo 
(Nanjing  Aeronautical  Institute) 

Abstract 

The  wall  lift  interference  parameters  on  ground  effects  for 
octagonal  closed  wind  tunnels  have  been  derived  using  image  vortex 
systems.  The  fillet  vortex  system  can  be  added  to  rectangular 
tunnel  vortex  system.  The  vortex  lattice  method  can  be  used  to 
determine  fillet  vortex  strength.  It  has  been  found  that  the 
wall  lift  interference  corrections  on  ground  effect  have  related 
to  not  only  the  wall  upwash  and  streamline  curvature  effects  but 
also  the  normal  gradient  of  the  upwash  velocity  at  the  horizontal 
tail . 

1.  Introduction 

The  interference  of  wall  lift  in  ground  effect  testing  has 

been  investigated  in  various  wind  tunnels  since  over  fifty  years 

ago.  The  image  vortex  systems  of  the  wing  model  in  the 

rectangular  wind  tunnel  are  two  groups  of  infinite  periodic 

vortex  series.  The  equation  for  the  wall  life  interference 

parameters  can  be  easily  derived  theoretically.  Batchelor  has 

designed  an  image  vortex  system  for  application  in  octagonal  wind 
[  1  ] 

tunnels 

He  has  based  on  the  rectangular  image  vortex  system  with  the 
addition  of  continuous  distribution  of  vortices  along  the 
rectangular  peripheries  to  determine  the  vortex  strength  using 


the  method  with  predetermined  pattern  and  to  calculate  the 
interference  parameter  of  the  wall  upwash  with  the  wing 
positioned  at  the  center  and  away  from  the  center  of  the  wind 
tunnel1'2'1. 

Fan  Jiechuan  has  applied  Batchelor's  method  to  calculate  the 

interference  parameter  of  wall  upwash  in  an  octagonal  wind 
[3] 

tunnel  .  He  has  also  extended  this  method  in  calculation  of 
the  ground  effects  and  obtained  several  equations  and  charts 
which  are  suitable  for  engineering  applications.  Garner  has 
indicated  in  reference  [4]  that  the  accuracy  of  the  calculation 
of  the  upwash  interference  parameter  using  Batchelor's  method  is 
considerably  poor  when  the  wingspan  of  the  wing  model  is  larger 
than  2/3  of  the  width  of  the  wind  tunnel.  Therefore,  we  have 
used  the  vortex  lattice  method  to  determine  the  fillet  vortex 
strength.  The  wall  upwash  interference  parameters  for  octagonal 
closed  wind  tunnels  calculated  by  this  method1^"1  are  quite 
consistent  with  the  accurate  solutions  obtained  .by  the  fixed  angle 
transformation  method1"^.  This  method  has  also  other  advantages 
such  as  it  is  simple  in  calculation  and  better  in  speed  of 
calculation. 

2.  The  Wall  Lift  Interference  Parameter  at  the  Wing 

A  testing  section  of  octagonal  tunnel  having  width  of  b, 
height  of  h,  and  fillet  length  of  t  is  considered.  The  wing  is 
assumed  to  be  located  at  the  center  of  the  tunnel  and  a  single 
horseshoe  vortex  having  strength,  r,  and  vortex  span,  1,  is 


representing  the  wing.  With  coordinate  system  oxyz  and  having 
the  x-axis  along  the  direction  of  incoming  flow,  assume  that  the 
attached  vortices  are  on  the  plan  of  x=0  and  the  free  vortices' 
are  extended  infinitely  to  the  downstream  and  parallel  to  the  x- 
axis.  In  the  octagonal  test,  every  fillet  is  divided  into  q 
elements  and  a  fillet  vortex  having  strength  is  placed  on 
every  1/4  point  ( z  j  ,  y  j )  .  The  3/4  point  (z^,  yNi)  t^le  contr°l 
point.  With  the  presence  of  wall  and  ground,  the  image  vortex 
systems,  in  respect  to  the  wall  and  the  ground  of  the  tunnel,  of 
the  horseshoe  vortex  of  the  wing  and  the  fillet  vortices  are 
shown  in  Figure  1. 

-  The  calculation  of  the  wall  lift  interference  parameter  for 
the  octagonal  wind  tunnel  can  be  carried  out  based  on  the 
calculation  for  the  rectangular  tunnel  having  the  same  width  to 
height  ratio 
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Figure  1.  The  Image  Vortex  Systems  of  the  Ground  Effect  Test 


with  the  additional  contribution  of  the  fillet  vortices.  With 
the  ground  effect,  the  image  vortex  systems  can  be  regarded  as 
four  groups  of  vertical  vortex  series  having  base  points  at  (z^ 
mb,  y j )  and  interspace  of  (h+2d),  where  ( z j ,  y j )  j  =  1,2, 3, 4  are 
the  coordinates  of  the  four  free  vortices  locating  immediately 

next  to  the  origin  of  the  coordinates,  and  m=0 ,  +1,  +2, . ,d 

are  the  heights  of  the  wing  from  the  ground.  Using  the  induced 
velocity  equation  for  the  vertical  vortex  series  having  equal 
space  and  same  direction  of  rotation,  the  equations  for 
calculating  the  average  wall  upwash  interference  parameter, 
at  the  wing  and  the  normal  derivative,  6yoR>  of  the  upwash 
interference  parameter  at  the  center  of  the  wing  in  the 
rectangular  wind  tunnel  can  be  derived  (see  reference  [5]  for 
detail ) . 

For  calculating  the  additional  effect  of  the  fillet  vortex, 
the  image  vortex  systems  of  the  wing  can  be  regarded  as  four 
equal  spaced  horizontal  vortex  series  locating  at  the  base 
points,  Ezj>  yj+n(h+2d)]  with  interspace  b,  where  n=0,  +1, 

+2, .  The  normal  velocity  at  any  control  point  on  the 

fillet  of  the  vortex  system  can  be  calculated  using  the  induced 
velocity  equation  for  the  horizontal  vortex  series  having  equal 
space  and  the  same  direction  of  rotation.  Similarly,  the  fillet 
vortex  systems  can  also  be  regarded  as  4  q  groups  of  horizontal 
vortex  series  based  at  Cz^,  yj+n(h+2d)]  having  interspace  b, 
equal  vortex  strength  and  the  same  rotation  direction.  (z^,  y j ) 
j  =  1,  2,...q  are  the  coordinates  of  the  vortices  on  tne  four 
fillets  locating  immediately  next  to  the  origin  of  the 


coordinates.  The  induced  normal  velocity  F.  .  at  the  control 

1  >  J 

point  of  the  vortex  series  having  unit  strength  can  be  obtained 
using  the  induced  velocity  equation  for  the  horizontal  vortex 
series.  With  the  condition  that  the  normal  velocities  on  the 
fillets  are  zero, 

±F,.IKI=C„  W-l.2-1) 

»-«  VI/ 

The  fillet  vortex  strength  can  be  obtained  by  solving  the 
above  group  of  linear  algebraic  equations.  The  upwash  velocities 
of  the  horizontal  clockwise  vortex  series,  having  base  points  at 
Czj,  yj+n(h+2d)]  and  interspace  b,  produced  near  the  wing  are: 


- ‘ilT  _ *  <2> 

. —  ch  ZL  Zy -*,-„<*+  2d)] -cos  *,) 

Using  the  above  equation  to  obtain  the  upwash  velocity  of  the 
wall  lift  interference  can  simplify  the  original  two-fold 
summation  formula  .down  to  single  summation.  This  will  give  high 
speed  and  accuracy  in  calculation.  Since  the  above  series 
converges  very  fast,  satisfactory  accuracy  can  be  obtained  with 
n=+3 . 

The  average  value  of  AV  .  can  be  obtained  by  integration  of 
equation  (2)  along  the  span  direction  on  the  wing  surface  (y=d). 
The  additional  average  upwash  interference  parameter  produced  on 
the  wing  by  all  fillet  image  vortex  systems  is: 


Obtaining  the  derivation  of  equation  (2)  with  respect  to  y,  tne 
normal  derivative  of  the  additional  induced  upwash  velocity  of 
the  fillet  systems  at  the  center  of  the  wing  (y=d,  z=0)  is 


bh‘  d(AK,,) 
2  n  £  dy 


(4) 


The  parametric  area  of  6^  and  6yo£  is  bh  and  the  reference  length 
of  the  latter  is  h.  The  normal  derivatives  of  the  average  wall 
upwash  interference  parameter  of  the  octagonal  wind  tunnel  and 
the  upwash  velocity  at  the  center  of  the  wing  are 

*-'*•+*>*  (5) 

Q 

^n=  ^  (6) 
respectively,  where  C  is  the  cross  area  of  the  test  section  of 
the  octagonal  tunnel.  The  detailed  equations  of  3  and  6yQand  the 
derivation  can  be  obtained  in  reference  [5], 

Figure  2  shows  the  dependence  of  the  wall  .interference 
parameter  6  on  the  span-to-width  ratio  1/b  and  the  relative 
height  above  ground,  d/h,  obtained  from  the  NH-2  low  velocity 
octagonal  wind  tunnel  of  the  Nanjing  Aeronautical  Institute  (the 
dimensions  of  the  octagonal  test  section  are  b=3M,  h=2.5M, 
t=0.4M).  It  indicates  that  the  wall  lift  interference  parameter  6 
decreases  rapidly  with  increasing  the  height  of  the  model  above 
i  ground.  The  value  of  6  with  ground  effect  is  smaller  than  that 

without  ground  effect.  The  difference  can  be  as  high  as  one 
order  of  magnitude. 


Figure  2.  The  Average  Wall  Upwash  Interference  Parameter  of 
the  NH-2  Wind  Tunnel 

1.  without  floor 

3.  1  he  Additional  Wall  Lift  Interference  Parameter  at 

Horizontal  Tail 

Since  the  horizontal  tail  is  located  downstream  of  the  wing, 
there  is  an  additional  upwash  velocity  increment  at  the  horizontal 
tail  compared  to  that  at  the  main  wing  under  the  influence  of  mapping 
fixed  vortex  and  free  vortex  systems  of  the  wing.  If  the  projected 
length  of  the  distance  between  1/4  chord  point  of  the  mean  aerodynamic 
chord  of  the  horizontal  tail  and  the  1/4  chord  point  of  M.A.C.  of  the 
wing  on  the  x-axis  is  L,  the  additional  upwash  velocity  produced  by 
the  wing  mapping  vortex  system  and  the  edge  mapping  vortex  systems  in 
the  octagonal  wind  tunnel  according  to  the  Biot-Savart  formula  is 
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where  o^,  a v  are  +1  or  -1  depending  on  the  direction  of  the 
vortex  o  o  =1  for  counter  clockwise  vortex  and  a  a  =-1  for  / 

4  V  MV 

clockwise  vortex  in  Figure  1,  =1,2...q  are  the 

coordinates  of  the  wing  and  the  fillet  vortices  respectively,  v 

and  y  are  integers  (Figure  1)  and  each  vortex  is  represented  by  a 

set  of  (m,v).  It  should  be  noticed  that  the  first  l  does  not 

m  ,  v 

include  (u,v)  =  (1,1),  (-1,1),  (1,-1)  and  (-1,-1).  They  are 
representing  the  horseshoe  vortex  and  the  images  with  respect  to 
the  ground  within  the  wind  tunnel.  Since  Equation  (7)  converges 
fast,  m=v=9  is  sufficient  for  calculation. 

From  Equation  (7),  the  additional  upwash  interference 
parameter  at  the  horizontal  tail  can  be  obtained  as^^ 


T<®«~  277  ,/r °)» 

and  the  normal  gradient  of  the  upwash  velocity  is^^ 
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(a)  interference  parameter 


(b)  interference  parameter 
1.  no  floor 

Figure  3.  The  Wall  Interference  Parameter  at  the  Horizontal 
Tail  in  the  NH-2  Wind  Tunnel  (l/b*0.5) 

Figure  3  shows  the  dependences  of  T6t  and  on  the 
relative  height  d/h  above  ground  and  the  relative  distance 
between  the  wing  and  the  horizontal  tail,  1/b,  in  the  NH-2  wind 
tunnel  with  a  span-to-width  ratio  of  l/b=0.5.  With  the  attack  ancle 
there  is  a  horizontal  distance  between  the  horizontal  tail  and  the 
center  of  the  wind  tunnel,  Ayt,  the  induced  increment  of  the 
attack  angle  at  the  horizontal  tail  is 


(10) 


where  is  the  uncorrected  wall  life  interference  parameter,  S 
is  the  area  of  the  wing.  Since  the  gradient  of  the  wall  upwash 
velocity  is  larger  at  the  horizontal  tail  and  the  position  of  the 


horizontal  tail  changes  with  the  attack  angle,  it  should  give 
attention  to  the  fact  that  the  errors  of  the  wall  lift  interference  corrections 
of  the  pitching  moment  and  longitudinal  static  stability  in  the 
ground  effect  testing  are  larger  if  the  second  term  in  Equation 
(9)  is  omitted^^. 

4.  Conclusions 

1.  The  wall  upwash  interference  parameter  with  ground 
effect  testing  is  smaller,  to  one  order  of  magnitude,  than  that 
without  ground  effect.  The  image  vortex  method  is  a  simple  and 
effective  method  of  approximation. 

2.  The  vortex  lattice  method  for  calculating  the  effect  of 

fillet  vortices  of  the  octagonal  wind  tunnel  is  simpler  and  more 
accurate  than  the  Batchelor's  predetermined  pattern  method.  /9 

3.  The  effect  of  the  normal  gradient  of  the  upwash  velocity 
at  the  horizontal  tail  should  be  considered  in  correcting  the 
wall  lift  interference  in  ground  effect  testing. 
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Abstract 

The  trail  lift  interference  parameters  on  ground  effects  for  octagonal 
closed  wind  tunnels  has  been  derived  using  image  vortex  systems.  The 
fillet  vortex  system  can  be  added  to  rectangular  tunnel  vortex  system. 
The  vortex  lattice  method  can  be  used  to  determine  fillet  vortex  strength. 
It  has  been  found  that  the  wall  lift  interference  corrections  on  ground 
effect  have  related  to  not  only  the  wall  upwash  and  streamline  curvature 
effeets,  but  also  the  normal  gradient  of  the  upwash  velocity  at  the  hori- 
aontal  tail. 
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